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cepts are  applied  to  generate  approximate  solutions  and 
filtering  algorithms*  and  the  well  known  extended  Kalman 
filter  eauations  and  higher  order  filtering  algorithms  are 
obtained  from  this  approach. 

The  concept  of  partitioning  the  measurements  is 
presented  and  shown  to  bring  advantages  in  computing  effi- 
ciency and  also*  for  nonlinear  measurements*  in  tracking 
accuracy.  A graphical  interpretation  of  the  action  of  Kalman 
filters  is  developed  and  provides  insight  into  the  impor- 
tance of  each  variable  in  the  filtering  process. 

Extensive  simulations*  designed  to  test  the  performance 
of  the  algorithms  developed*  are  presented  in  graphical  form 
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SYMBOLS  AND  CONVENTIONS  FOR  TEXT 

1 - Underlined  small  letters  represent  vectors#  thus# 
x.  denotes  the  vector  x. 

2 - Capital  letters  are  used  to  represent  matrices. 

3 - A caret  indicates  an  estimated  value#  e.g.#  x 
denotes  an  estimate  of  x. 

- Time  points  are  represented  by  t^  or#  when 
between  parentheses#  simply  by  k»  e.g.#  x_(  k ) denotes  the 

value  of  x at  time  t,  . 

— k 

5 - Estimates  specify  the  time  of  estimation  and  the 
amount  of  information  used#  e.g.#  x(k,'j)  denotes  the  esti- 
mate of  _x  at  time  t^takinq  into  account  all  the  measurements 
uo  to  and  including  t^  . When  k = j only  one  time  point  is 
indicated#  thus#  x^(k)  denotes  the  estimate  of  at  time  tf 
taking  into  account  all  measurements  uo  to  and  including  t^ . 

6 - Probability  density  functions  are  represented  by 
the  small  letter  p#  i.e.  Px|y(*'v)  denotes  the  conditional 
density  of  x given  y.  when  no  confusion  is  oossible#  this  is 
simplified  to  p(x!y). 

7 - The  expectation  operator  is  represented  by  E f 1; 
the  variance  by  Vart  1. 
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at 


a(k) 
o 8(k) 

“b(k) 
o f(  k) 

b ( k ) 

bi(  k) 

F(k) 


expect'  value  of  the  estimated  state  vector  at 


random  rate  of  change  of  target  speed. 


random  rate  of  change  of  target  heading. 


random  rate  of  change  of  frequency  emitted. 


heading  of  target  at  time  t^. 


bearing  measurement  by  buoy  i 


second  moment  of  the  state  estimate  around  expected 


value 


d iC  k ) 


kj 

e ( k ) 
e ( k ) 


f (k) 
o 


f Ck) 


G(k) 
H ( k ) 


PCk ) 


a(k). 

average  speed  of  sound, 
distance  of  target  to  buoy  i. 

Kronecker  delta;  - 1 if  k - jf  - 0 if  k 
estimation  error  vector. 

exoected  value  of  the  estimation  error, 
rest  frequency  emitted  by  the  target, 
frequency  measured  by  buoy  i. 
gain  matrix, 
observation  matrix, 
compression /expansion  factor 
estimation  error  covariance  matrix. 


= i • 


Q ( k ) state  excitation  covariance  matrix. 


R(k)  measurement  noise  cavariance  matrix. 


s ( k ) speed  of  target. 


s^ik)  relative  velocity  of  target  towards  buoy  i. 


a (k)  standard  deviation  of o (k) 
s s 


o,(k)  standard  deviation  of  (k) 
b b 


a^(k)  standard  deviation  of  (k). 


time  delay. 


T time  delay  difference  between  the  signals  received 

by  buoys  i and  j . 


time  spent  in  prediction. 


P 

T time  scent  in  estimation, 

e 

v(k)  vector  of  q random  measurement  noise  signals. 

jw(k)  vector  of  m random  forcinq  inputs. 

_x(k)  state  vector  of  dimension  n. 

x(k)  x-position  of  target. 

x^(k)  x-oosition  of  buoy  i. 

x_C  k ) estimated  state  vector. 

y/ k ) tn  Cn  + 3 ) /2] -di  mens  i ona  1 vector  containing  the  state 
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variables  and  the  distinct  components  of  the  matrix  P< 


v ( k ) 

y-pos i t i on 

of  target. 

Vi  (k) 

y-pos i t i on 

o f buoy  i . 

z(k> 

vector  of 

q measurements. 

Zk 

set  of  all 

measurements  up  to  and  including  time  t,  > 

k 

i .e. 

tz(0),  z(t). 

. . . f Z ( k ) 1 . 
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I.  INTRODUCTION 


This  work  was  oriented  to  the  problem  of  optimally 
estimating  characteristics  and/or  parameters  of  a certain 
class  of  nonlinear/  dynamic/  stochastic  systems  from 
observed  output  sequences  which  are  noise  corrupted  non- 
linear functions  of  the  system  states. 

In  ©articular/  attention  was  directed  to  the  problem  of 
real-time#  short-range  optimal  localization  and  tracking  of 
a submarine  target  by  passive  acoustical  means.  The  most 
accurate  and  reliaole  estimates  of  the  target's  parameters 
(position/  heading/  soeed)  are  sought  by  optimally  process- 
ing the  measurements  obtained  from  the  acoustic  signals 
transmitted  by  the  target  itself  and  received  by  special 
sonobuoy  s . 

In  the  unclassified  literature  there  are  presently 
available  methods  of  passive  tracking  of  submarine  targets. 
Reference  UJ  utilizes  dooo 1 er-sh i f t ea  frequency  measure- 
ments obtained  from  a group  of  sensors;  [21  uses  bearing  and 
dopp 1 e r-sh i f t ed  f reauency  measurements  obtained  from  one  or 
more  sonobuoys.  These  methods  provide  very  good  first 
approx i mat i ons  to  the  solution  of  the  problem  and  many  of 
their  concepts  and  notation  are  used  in  this  work. 


As  a first  step*  an  attempt  was  made  to  produce  a 
comprehensive  study  of  all  the  information  available  in  the 
acoustic  signals  picked  up  by  the  sonobuoys  and  its  useful- 
ness in  the  estimation  process. 

The  estimation  algorithms  were  developed  taking  into 
account  the  advantages  of  processing  information  as  soon  as 
it  becomes  available.  Flexibility  to  move#  remove  and 
include  sensors  ana  measurements  at  any  time  was  obtained. 
Constraints  were  included  to  account  for  limitations  of 
practical#  inexpensive  and  presently  available  sonobuoys  and 
computing  equipment. 

In  Chapter  II  the  problem  is  described  in  mathematical 
terms  using  state  space  techniques  and  the  initial  assump- 
tions and  constraints  are  given.  The  model  used  for  the 
target  is  similar  to  the  model  presented  in  [21#  but  with  a 
different  set  of  states.  The  measurements  obtained  from  the 
signals  received  by  the  sonobuoys  are  characterized  and 
their  functional  relationships  to  the  parameters  of  the  tar- 
get are  ana  1 y zed . 

Chapter  III  discusses  the  general  theoretical  solution 
for  the  problem  and  the  difficulties  of  implementing  it. 
Approximate  solutions  are  then  sought  ana  processing  equa- 
tions are  developed.  The  special  vectors#  matrices  and 
relations  characteristic  of  the  tracking  problem  are 
prepared  for  the  application  of  the  filtering  eauations. 


la 


In  Chapter  IV  the  concept  of  partitioned  measurements 
is  introduced  and  analyzed.  The  advantages  in  accuracy*  com- 
puting efficiency  and  filter  flexibility  are  stressed. 

Practical  and  graphical  interpretations  of  the  estima- 
tion process  help  in  visualizing  nonlinear  problems  and 
motivate  the  adoption  of  iterative  techniques.  Increased 
accuracy  and  possible  convergence  improvements  are  shown  to 
result  from  this  approach. 

Chapter  V introduces  the  interactive  computer  program 
that  allowed  the  application  of  the  ideas  describee  in  the 
previous  chapters  to  the  tracking  problem.  Selected  simula- 


tions  are 

di scussed 

and  presented 

in  graphical 

form. 

This 

research 

covers  the 

subject 

i n 

a muc  h 

more 

comprehensive  way 

than  earl i er 

works* 

m * 

121*  131 

. The 

results  obtained  by  using  the  methods  and  ideas  collected 
and  developed  in  this  study  show  many  ways  of  improving  the 
tracking  accuracy*  of  increasing  computing  efficiency*  of 
reducing  divergence  problems  and  of  obtaining  faster  conver- 
gence . 


The  importance  of  having  the  buovs  placed  in  adeauate 
positions*  of  havinq  an  adequate  model  for  the  target*  and 
of  having  freouent*  varied  and  precise  measurements  is 
st  ressed. 
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II.  PROBLEM  FORMULATION 


V 

A.  THE  ESTIMATION  PROBLEM 

A general  nonlinear#  stochastic#  dynamic  system  observed 

* 

by  a set  of  nonlinear,  nonstationary  measurements  can  be 
represented,  for  our  purposes,  in  state-space  discrete  form 
by  the  equations 

x(k+l)  = f(x(k),w(k),t.  ,t,  ) (2.1) 

— — — — k k+1 

z_(k)=_h(x_(k),v_(k)»t,  ) (2.2) 

where  x^(k)  is  the  state  vector  of  dimension  n, 

v,(k)  is  the  vector  of  the  m random  forcing  inputs, 
^(k)  is  the  measurement  vector  of  dimension  a, 
v^(k)  is  the  vector  of  p random  measurement  noises, 

_f(  ) and  h(  ) are  general  vector  functions. 

The  general  estimation  problem  consists  of 
---  knowing  the  statistical  characteristics  of  w(k)  and 
jv(k)  and  having  a statistical  representation  of  the  initial 
conditions  x^(O)  of  the  system, 

— processing  the  observations  2^°^'  z.(l ) , ...,  £(j) 

in  a statistically  optimal  way  to  obtain  the  best  estimate, 
in  some  sense,  of  the  system  state  x(k)  at  time  t,  . 

If  k > j the  processing  is  called  prediction,  if  k < j 
it  is  called  smoothing  and  if  k = j it  is  called  filtering. 


In  this  work  some  departures  from  the  completely  general 
nonlinear  case  were  taken.  The  first  assumption  is  that  the 
measurement  noise  is  additive  and  that  the  system  equations 
are  linear  with  respect  to  the  random  input  vector  w_.  Equa- 
tions (2.1)  and  (2.2)  can  then  be  written  as 


x(k  + l) 


f ( x ( k ) » t,  » t ) ♦ g(x(k)*t  *t  ) . w ( k ) 

k k+1  - - k k+1  - 

(2.3) 


z^k  + 1)  s h_(j<(k)»t^)  + v^( k ) 


(2. a) 


where  f_(  ) and  M ) represent  general  vector  functions  and 
g(  ) a general  matrix  of  functions. 


Secondly*  the  measurement  noise  and  random  forcing  input 
are  assumed  to  be  uncorrelated  in  time*  zero  mean*  discrete 
Gaussian  sequences*  independent  of  one  another  and  indepen- 
dent of  the  initial  condition  of  the  state  of  the  system, 
i .e . * 

E (w(k)l  = 0 Etv(k)l  = 0 (2.5) 

E(w(k)wT(j)l  = Q ( k ).  6,.  Elv(k)vT(j)l  = R(k).  6,. 

- - kj  - - kj 

E lw(k)v_T(  j ))  = 0 EU(0)wT(k)l  = 0 E l v ( k )^T  ( 0 ) 1 = 0 


1 i f k.  = j 

0 i f k i ) 
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B.  THE  PASSIVE  TRACKING  PROBLEM 


Many  practical  situations  can  be  formulated  as  estima- 
tion problems  and  characterized  by  the  equations  described 
above.  The  situation  that  motivated  this  study  assumes  a 
submarine  target  following  (most  of  the  time)  an  approxi- 
mately constant  speed  and  constant  heading  path  in  a field 
of  passive  sonobuoys. 

The  target  un i n t en t i ona 1 1 y emits  acoustic  signals  which 
are  picked  up  by  hydrophones.  The  transduced  signals  are 
sent  up  by  wire  to  the  buoys  and  then  f reguency-modu 1 ate  VHF 
carriers  which  are  transmitted  by  the  buoys  and  receiveo  at 
a nearby  ship  or  aircraft.  After  the  recovery  of  the  signal# 
data  processing  equipment  generates  measurement  values  which 
contain  numerical  information  about  the  target  parameters. 

The  first  stoo  in  the  analysis  of  the  signals  collected 
by  the  sonobuoys  is  the  attempt  to  detect  the  presence  of  a 
real  target. 

After  a tarqet  is  detected#  frequency  spectrum  informa- 
tion is  added  to  intelligence  data  in  an  attempt  to  classify 
« 

it.  Other  measurements  and  information  from  any  other  source 
are  also  used  to  obtain  a first  approximation  to  the 
target's  parameters  (position#  heading#  speed)#  generally 
through  the  use  of  least  mean  square  techniques  (1)#(3). 

The,., third  phase  of  the  process#  if  a wartime  condition 
doesn't  exist#  is  the  tracking  of  the  target.  That  is#  the 


determination  of  its  oath  while  reducing,  if  possible,  the 
initial  uncertainties  about  its  parameters,  in  real-time, 
through  the  optimal  processing  of  the  sparse  information 
provided  by  the  sonobuoys.  The  passive  characteristic  causes 
severe  limitations  but  is  necessary  to  avoid  revealing  to 
the  target  the  fact  that  it  is  being  observed. 

This  investigation  was  devoted  to  the  tracking  phase  and 
it  was  assumed  that  the  main  characteristics  of  the  buoys  as 
well  as  their  positions  are  known. 

C.  THE  MODEL  OF  THE  TARGET 

The  problem  is,  o*  course,  three-dimensional.  Neverthe- 
less the  measurements  are  not  a direct  function  of  the 
target's  deoth  but  of  the  difference  in  depth  between  the 
target  and  the  hydrophones.  Anticipating  the  accuracy  of  the 
measurements  and  the  precision  of  the  data  processing  eouip- 
ment  and  algorithms,  it  is  justifiable  to  consider  only  two 
dimensions  initially,  adding  depth  as  the  third  dimension 
whenever  necessary  and  with  no  conceptual  difficulty. 

The  target's  path  is  freauently  one  with  nearly  constant 
SDeed  and  heading,  disturbed  by  currents  and  random  thrust 
and  control  variations.  Intentional  maneuvers  which  have  no 
evasive  purposes  are  normally  simole  and  smooth. 

A basic  plant  having  as  states  two  position  coordinates 
(x  and  y),  heading  (b)  and  speed  (s)  of  the  target,  describ- 


i ng  a constant  velocity  path, 


as  shown  in  Fig  1 , 


seems 


- - ■ --  | -™'~'  H ■ 


adequate . 

Random  forcing  inputs  should  be  considered  to  account 
for  noise  processes  and  imperfections  of  the  model*  such  as 
target  maneuvers. 

As  suggested  by  12)  these  random  forcing  functions  can 
be  approximately  represented  by  two  independent  zero-mean* 
pi ecew i se-const ant  random  rates  of  change*  a and  o * act- 

S D 

ing  on  the  speed  and  heading  of  the  target. 

The  formulation  of  this  basic  plant  in  discrete  state- 
space  form*  similar  to  Eauation  (2.3)*  can  be  obtained  in 
the  following  way 


s ( k + 1 ) - s(k)  ~ 
b ( k + 1 ) - b ( k ) = 
x ( k+ 1 ) - x(k)  = 

ft. 


k+1 


a (t)  dt  = (t  - t ) . a C k J 
s k+1  k s 


ft 


k+1 


uk 

fCk+l 


a ( t ) dt  = (t,  , - t ) .a  (k) 

b k+1  k b 


s(t)  cos  bCt)  dt  = 


k+1 


(s(k)  + Ct  - t,  )a  ( k ) 1 coslb(k)  + 
k s 


* . 

t * 


♦ (t  - tfc)ab(k)l  dt  = 


k+1 


s ( k ) cos  b(k)  dt  ♦ 


k+1 


Ck 

ck+l, 


(t  - t,  )a  (k)  cos  b(k)  dt  - 
k s 


'( t - t,  )a.(k)  s(k)  sin  b(k)  dt  ♦ ... 
k b 


ck 


s (t.  .1  - t,  ) s(k)  cos  bCk)  t 

k+1  k 


21 


J 


♦ 


J. 

2 


(t 


k+1 


( a ( k ) cos  b ( k ) - 
8 


- a (k)  s ( k ) sin  b(k)]  ♦ HOT 
b 


where  HOT  represents  Higher-Order  Terms. 


In  the  same  wayf 


(2.6) 


y (k  + 1 ) 


- y ( k ) 


(t 


k+1 


s(t ) sin  b(t ) dt  = 


(t.,-  - t)  s(k)  sin  b(k)  + 

k+1  k 

♦ (t  - t )2  (a  (k)  sin  b(k) 

2 k+1  k s 

♦ a (k)  s(k)  cos  b(k)J  + HOT 

D 


♦ 


: 


(2.7) 


If  ltk+l  ’ tk)as(k)  and  (tk+l  ' *k  )ab  U)  ar*  SUfH" 
ciently  small  so  that  the  higher  order  terms  in  Equations 

(2.6)  and  (2.7)  can  be  neglected#  one  has#  in  vector  form 


x_(k  + l)  = # t.  . ) ♦ g_(x(k)#t  #t  ).w(k) 

iCt  X k k+1 


where 


x(k)  = 


x ( k ) 
y ( x ) 
s ( k ) 
b ( k ) 


(2.8) 


a (k) 
s 

a ( k ) 


w ( k ) 


(2.9) 


f ( x ( k ) , t , # t ) 

~ - k+1  k 


map  mm 


(2.10) 


x ( k ) ♦ (t 

- t ) 
k+1  k 

s(k) 

cos 

b ( k ) 

y ( k ) ♦ ( t 

- t ) 
k+1  k 

s ( k ) 

s i n 

b ( k ) 

s ( k ) 
b(k) 

and 

g(x(k)#t#t)=  (2.11) 

— — k+1  k 


2 ( * k+1 

2 

- tfc)  cos 

hlk>  - i'Vi 

- 

t,  )2  s(k) 
k 

sin  b ( k ) 

I(t 

2 1 1 k+1 

2 

* tk)  sin 

b(k)  i'Vi 

- 

t,  )2  s ( k ) 
k 

cos  b ( k ) 

(tk+l  " *k 

) 

0 

- 

0 

(t 

k+1  “ *kJ 

when 

f r equenc  y 

measurements  are 

to 

be 

processed 

di rect 1 y » 

as  shown 

i n Sec  tion  1 1 » D # 1 » o , 

an 

addi t i onal  state 

variable  must  be  included  ---  the  rest  f reauency  f emitted 

o 

by  the  target. 


This  frequency  is  assumed  to  be  aporo x i ma t e 1 y constant 
with  a small  random  disturbance  a £ that  is  assumed  to  be 
zero-mean#  pi ecewi se-constant » and  independent  of  ag  and  a b . 

The  expanded  state  variable  formulation  is  given  by 


x e ( k ) = 


x ( k ) 


f (k) 
o 


(2.12) 


2« 


we(k)  = 


(2.13) 


_xe(k  + l)  = i (£.  (k),tk+1  ,tk)  ♦ 


+ g e(  xe  ( k ) > t , , » t . ) . we(  k ) 

— — k+1  k — 


(2.1*0 


*here 


le(*_e(k),tk+1,tk)  * 


-(- C k 1 ' *k+l  ' lk  ) 


f (k) 
o 


( 2 • 1 5 ) 


ge  ( x e ( k ) , t ft  ) s 

— — k+1  k 


3(x(k)ftk+1  ,tk  )j 
1 

0 

1 

1 

1 

0 

— 

o o ! (t 

- t , ; 

i 

k+l  k 

(2.16) 

D.  MEASUREMENTS 

The  sonobuoys  can  perform  one  or  both  of  the  following 
tasks : 

---  pick  up  underwater  sound  signals 

---  indicate  the  approximate  direction  of  the  vectorial 
sum  of  the  received  acoustic  signals. 


m 


The  target  emits  a signal  r(t)  which  is  distorted  and 
modulated  by  the  propagation  medium  between  the  source  and 
the  hydrophone*  and  by  extraneous  sound  sources  present  in 
the  ocean. 

When  all  of  these  disturbances  are  of  relatively  mild 
strength*  the  signal  received  by  a stationary  buoy  i will  be 
approximately  an  attenuated  and  delayed  version  of  the  pre- 
viously emitted  signal,  i.e.* 


r (t)  = a (t).r(t  - t (t)) 
i i i 


(2.17) 


where  a.(t)  is  an  attenuation  factor  and  t ^ ( t ) is  the  time 
delay  given  by 


r (t)  = d (t  - T (t))  / c 
i i i 


(2.18) 


that  is*  the  distance  d^  between  the  target  and  the  buoy  at 
the  time  of  emission  divided  oy  the  average  velocity  of 
sound  in  the  medium. 

Since  the  target  velocity  is  much  smaller  than  the  sound 
velocity*  for  small  distances  and  delays  one  has 


d (t  - x (t))  = d (t)  + s (t).T  (t) 
i i i i i 


(2. 1*1) 


where  s.(t)  is  the  relative  speed  of  the  target  toward  buoy 
i as  shown  in  Fig  2. 

Thus,  from  Eguations  (2.18)  and  (2.19)  one  has 


t i ( t ) = di(t)  / (c  - s (t)> 


(2.20) 


2b 


I 


and  the  received  signal  is  qiven  by 


r i ( t ) = ai(t).rlt  - di(t)/(c  - s A ( t i ) J 


(2. 21) 


If  this  signal  is  recorded  from  t = tQ  to  t * tQ  + T, 
one  has  recorded 


r It  + t ' ) = a (t  tt’l.rlt  ♦ t ' - 
i o i o o 


-d(t  * t')/(c  • s (t  ♦t'))J 
i o i o 


° < t ' < T 

If  during  this  oeriod  T the  relative  target  speed 
towards  buoy  i doesn't  change  significantly# 


s ( t +t')=s(t) 
i o 10 


d(t  +t')=d(t)-t’.sCt) 
i o i o i o 


and  hence 


*(t  + t ' ) = a C t ♦ t ' ) . r ft  -d(t)/(c-s(t)J* 

i o i o '—o  i o l o 


+ t'.ll  + s (t  )/(c  - s (t  ) ) f) 

i o i o 


r (t  + t')  = a (t  ♦ t'J.rlt  - a (t  )/Cc  - s (t  ))  ♦ 
io  io  oio  io 


* c.t'/(c  - s (t  ))) 
i o 


(2.22) 


In  equation  (2.22),  the  second  term  inside  the  brackets 


i s a 


fixed  time  delay.  The  third  term  is  a variable  time 


delay  and  can  be  seen  as  a comoress i on/expans i on  term,  show- 


ing  the  different  value  of  the  variable  time  for  the  origi- 
nal and  the  received  signals. 

If  the  emitted  signal  had  a single  freguency  component  wQ 
one  would  have 

r (t  ♦ t ' ) = a (t  ♦ t').cos(w  (t  ♦ t')) 
i • o io  i o 

where 

cos(w  (t  + t ' ) ) = cos  ( it  ♦ w .t')  • 
i o yio  i 

• 

= costw  ,t  - w ,d  (t  )/(c  - s (t  ))  ♦ 
oo  o i o i o 

♦ w .c.t'/Cc  - s (t  ))1  = 
o i o 

= cos  t <p  ♦ w .c.t'/Cc  - s (t  ))] 
io  o i o 

and  the  received  frequency  would  be 

w = w.c/(c-s) 
i o i 

fsf.c/Cc-s)  (2.2 3 ) 

i o i 

The  signals  collected  by  various  sonobuoys  are  sent  to  a 
ship  or  aircraft  normally  eauipped  with  equipment  and  skills 
to  extract  the  information  described  below. 

I . Isolated  Buoys 

By  processing  the  signals  collected  by  each  buoy 
alone»  one  can  obtain  the  following  noisy  measurements: 


a.  Bearing  Indication 

Bearing  is  given  directly  bv  the  buoys  or 
preprocessed  (filtered)  to  reduce  noise  effects. 

Since  the  bearing  is  obtained  from  a vectorial 
suit*  it  can  be  reasonably  accurate  only  when  the  acoustic 
noise  is  approximately  omnidirectional  or  much  weaker  than 
the  target  signal.  For  some  sensors  this  vectorial  sum  can 
be  limited  to  a selectable  frequency  band*  thereby  easing 
the  noise  problem. 

The  re  1 at i onsh i ps  between  a bearing  indication 
obtained  by  a buoy  i and  the  states  of  the  target  is  shown 
in  Fig  3 and  given  by 

b^k)  = arctanQy(k)  - y.(k)l/(x(k)  - x^(k)Q  + v(k) 

(2.24) 


where  v(k)  must  account  for  all  the  noise  aisturoinq  this 
measurement.  including  acoustic  noise,  transducer  inaccura- 
cies. preprocessing  noise  and  errors  due  to  the  finite  speed 
of  the  sound  in  the  water.  This  last  factor  is  caused  by  the 


time  delay  between  the  emission  of  a sound  by  the  target  and 
reception  at  the  hydrophone.  At  the  end  of  this  period  the 
target  and  the  hydrophone  have  slightly  different  relative 
positions. 

According  to  (2)  tyoical  accuracies  of  _t  5 
degrees  are  common  for  strong  signals  and  inexpensive  DIFAR 
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or  less  can  he  obtained  with  arrays  of  hydrophones  and 
proper  preprocessing. 

b.  Frequency  Indication 

By  analyzing  records  of  length  T seconds  of  the 
acoustic  signals  picked  up  by  the  sonobuoys  with  Fast 
Fourier  Transform  algorithms*  one  can  obtain  approximate 
frequency  spectrum  descriptions  with  a resolution  of  1/T 
Hertz. 

By  detecting*  measuring  and  possibly  tracking 
[4]  some  of  tne  strongest  frequencies  contained  in  the  sig- 
nal* very  useful  data  about  the  target  parameters  is 
obtained  since  the  received  frequencies  are  functions  of  the 
originally  emitted  frequencies*  the  speed*  the  heading  and 
the  position  of  the  target  * as  shown  in  the  relationship 

f (k)  = f (k).c  / Cc  - s (k))  + v(k)  (2.25) 

i o i 

where 

s^k)  = -s  ( k ) .cos  lb  ( k ) - b^  ( k ) 1 (2.26) 

is  the  relative  velocity  toward  buoy  i*  as-shown  in  Fig  3. 

The  term  v(k)  now  has  to  account  also  for  the 
errors  introduced  in  the  various  paths  of  the  signal  from 
the  target  to  the  ship  or  aircraft*  and  for  variations  in 
t a rget -sensor  geometries  during  the  record  time  T. 
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As  suggested  in  (11*12)  and  13) * inaccuracies 

ranging  between  0.5  and  0.01  Hertz  can  be  obtained  in  most 

situations  for  values  of  f of  the  order  of  hundreds  of 

o 

Hert  z . 


c.  Signal  Strength  Indication 


The  strength  of  the  signals  varies  not  only  with 
the  strength  of  the  originally  emitted  sound  but  also  with 
the  distance  between  the  target  and  the  buoy.  The  signal 
• strength  also  depends  on  the  random  aspects  and  the  direc- 
tivities of  the  target*  the  hydrophones*  and  of  the 
transmitting  and  receiving  VHF  antennas.  Random  absorption 
and  scattering  in  the  various  media  from  the  target  to  the 
processing  equipments  may  also  have  a significant  influence 
on  the  signal  presented  to  the  processors. 


A mathematical  formulation  of  all  these  effects 
would  be  extremely  difficult  task  and  would  obscure  the 
desired  distance  information.  For  this  reason  the  informa- 
tion contained  in  the  strength  of  the  collected  signals* 
which  can  be  used  by  an  experienced  and  skilled  audio  opera- 
tor* was  not  incorporated  in  this  study. 


Attenuation  factors  are  thrown  away  by  setting 
the  strength  of  the  received  signals  at  a good  working 
level  . 
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Two  op  More  Buoys 


Considering  the  signals  received  by  two  or  more 
buoys*  indications  can  be  obtained  which  may  be  used  in 
addition  to*  or  instead  of*  the  measurements  available  from 
a single  buoy. 

a.  Frequency  Difference  Indication 

To  process  frequency  measurements  obtained  from 
one  buoy  requires  inclusion  of  the  rest  frequency  as  one  of 
the  states.  This  is  the  case  because  this  measurement  is 
very  sensitive  to  errors  in  the  assumed  rest  frequency* 
i .e.  # 

f . = f .c  / (c  - s . ) 
i o i 

A f . - Af  .c  / (c  - s.)  = Af 
io  i o 

Having  the  frequencies  received  by  two  buoys  and 
taking  their  difference  one  gets 

f . " f . = f .c/(c  - s ) - f .c/Cc  - s.)  = 
i j o i o j 

= fQ.c.(si  - s^)  / C ( c - s ^) . (c  - S^)l  = 

= f o. (s  ± - s )/c  (2.27) 

The  sensitivity  of  this  difference  to  errors  in 
the  assumption  of  the  rest  frequency  is  given  bv 

A(f  - f ) = Af  .(s.  - s )/c 
i j o 1 j 

3U 

J 


op  is  reduced  by  a very  large  factor  since  c is  of  the  order 
of  1500  m/s  and  (s^  - s^)  is  normally  very  small. 

This  way  * if  one  has  a reasonable  approx i mat i on 
for  f t instead  of  processing  two  measurements  in  a five- 
state  plant  one  can  process  only  one  f reguenc y-di f f erence 
measurement  in  a four-state  olant  with  greatly  reduced  com- 
puting time  and  hopefully  not  a detectable  reduction  in 
accuracy . 

The  variance  of  the  errors  in  each  measurement 
must  be  added  to  give  the  variance  of  the  error  in  the  fre- 
quency difference  measurement. 

b.  Frequency  Ratio  Indication 

If  one  divides  the  frequencies  measured  at  two 
sonoouoysr  a new  relationship  is  obtained  that  retains  the 
information  on  oosition*  heaainq  and  speed  of  the  target* 
but. is  independent  of  the  rest  frequency  emitted 

f,  / f = tf  . c / ( c - s,)).t(c  - s.)/(f  .c)J  = 
i j o i jo 

= (c  - s ) / (c  - s ) (2.26) 

i i 

or 

f ^ / f , : 1 ♦ (s^  ■ Sj)/c  . tl  ♦ si/c  ♦ (s^/c)2  + ...)  = 

s 1 ♦ (st  - s ) / c (2.29) 
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As  with  the  frequency  difference  indication# 
this  measurement  can  be  used  in  a four-state  plant  instead 
of  using  two  frequency  measurements  in  a five-state  plant. 

If  this  ratio  is  generated  by  the  division  of 
the  frequencies  obtained  from  the  analysis  of  each  of  the 
signals#  the  characteristics  of  the  frequency-ratio  measure- 
ment noise  is  complicated  and  state  dependent.  If  however 
this  relation  is  obtained  as  described  in  the  following  dis- 
cussion# this  uncomfortable  situation  can  be  avoided. 

c.  Time  Delay  Indication 


As  shown  at  the  beginning  of  this  section  on 
measurements#  the  signal  collected  by  a sonobuoy  is  a noisy# 
delayed#  attenuated  and  comoressed/exoanded  reproduction  of 
the  original  signal  emitted  by  the  target.  The  delay  is  due 
to  the  finite  time  the  signal  takes  to  reach  the  buoy  and 
the  compression/expansion  is  caused  by  the  variation  of  this 
delay  with  time  (doooler  effect). 


The  signals  received  by  two  sonoDuoys  have  dif- 
ferent noise  contributions  and  different  delay  and 
compression/expansion  factors  even  after  their  strengths  are 
eoua 1 i zed. 


IS: 

i 

* ^ 

M 


Supoose  one  recorded 
buoys  i and  j from  t 
related  to  the  original  signal  as 
i . e . # 


t to  t 

0 
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the  signals  received  oy 

= t t T.  The  signals  are 

o 

shown  in  Equation  (d.22), 


•/ 


r . ( t +•  t * ) = a.(t  ♦ t ' ) . r 1 1 - ci  (t  )/(c  - s (t  ))  ♦ 

i o io  oio  io 


♦ c . t * / (c  - s ( t ))1 
i o 


(2.30) 


r .( t ♦ t * ) = a ( t t t'J.rlt  - d.(t  )/(c  - s (t  ))  ♦ 
j o jo  ojo  jo 


♦ c . t ' / (c  - s ( t ))] 
j o 


(2.31  ) 


If  one  now  amplifies  these  two  signals  to  a good 
working  level » correlate  them  by  evaluating 


p ( t)  = 


fT 

= Pi 


(t  ♦ t ' ) . r ( t ♦ t ' + t l.dt' 
o jo 


(2.32) 


and  look  for  the  value  of  t that  maximizes  p (t)  one  finds 
that  : 


(1)  if  s (t  ) : s (t  ) and/  conseouent 1 y / the  fre- 
10  jo 

quency  shifts  are  aoout  the  same  in  both  signals/  then  r (t 

i o 

+ t * ) will  be  sirndy  an  aooroximate  delayed  or  advanced 

replica  of  r (t  t t')/  as  can  be  seen  from  the  above  eaua- 
3 o 


tions. 


In  this  case  the  maximum  of  p ( x ) is  easy  to 


find  and  corresponds  to  the  value  x at  which 

max 


t - d . ( t )/(c  - s ft  ))  ♦ c . ( t 1 + T ) / ( c - s (t  )) 
oio  i o max  i o 

= t - d (t  )/(c  - s (t  ))  ♦ c.t'/(c  - s (t  )) 
ojo  jo  jo 


solving  for  x yields 
max 


x = (d  (t  ) - d (t  )]  / c 
max  1 o jo 


(2.3.J) 
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(2)  if  s.(t  ) 1 s.(t  )»  o ( t ) does  not  have  a 
io  jo  K 

simply-determined  maximum  value  because  the  two  signals  have 
frequency  components  whose  phases  are  shifted  by  different 
amount  s . 

If  ones  now  change  the  time  scale  of  one  of 

the  signals*  say  r.(t  + t')»  one  has 

j o 


r'(t')  = r. (t  ♦ m . t ' ) = 

j jo 


= a.  .rflt  - d.Ct  )/(c  - s,(t  1)1  + 
j L o jo  jo 


♦ C .m . t ' / (c 


(2.34) 


Comparing  this  equation  with  Equation  (2.30) 
one  sees  that  if  m can  be  found  such  that 


m.c/(c  - s.(t  ))  = c/(c  - s (t  )) 
jo  i o 


or 


m = Ic  - s.(t  )1  / [c  - s.(t  )1 
jo  i o 


(2.35) 


one  again  has  case  (1)  and  the  correlation  function  will 
give  a maximum  at  approximately 


t = Id . ( t ) - d ,(t  ))  / c (2.36) 

max  io  jo 


It  is  interesting  to  note  that  Equation  (2.35) 
is  the  same  as  Equation  (2.20).  The  compress i on/exoans i on 
factor  is  thus  the  ratio  between  the  frequencies  received  Dy 
buoys  i and  j . 
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The  p (-wrr)  function#  if  the  two  recorded  sig- 
nals were  broadband#  of  long  duration  and  with  high  S/N 
ratios#  would  appear  as  shown  in  Fig  U. 

Low  signal  to  noise  ratios  could  hide  the  real 
maximum  and  create  false  ones#  narrowband  or  coherent  sig- 
nals would  present  ambiguities  with  the  creation  of  many 
close  maxima. 

By  considering  the  availability  of  equipments 
and  techniques  to  obtain  the  prooer  elimination  of  the 
difference  in  frequency  shifts  in  the  recorded  signals  # 
this  study  assumed  the  possibility  of  having  time  delay 
measurements  available  to  provide  information  about  the  tar- 
get position  accordina  to  the  relationship 

T.  (k)  = Id  (k)  - d.Ck)J  / c + v(k)  U.  37) 
ij  i j 


where  v(k)  must  account  for  the  errors  caused  by  the  dif- 
ferent noise  contributions  in  each  signal,  for  the  accuracy 
and  orecision  of  the  cross-correlation  algorithm  or  device# 
and  for  variations  in  target-sensors  aeometry  during  the 
recording  of  the  signals. 

Inaccuracies  anywhere  from  0.5  sec  down  to  a few 
milliseconds  seem  very  reasonable  to  expect  in  many  cases. 


d.  Siqnal  Strength  Difference 

For  the  same  reasons  already  explained  in 
II#D#l#c  the  information  contained  in  the  difference  of 
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III.  THEORETICAL  SOLUTION 


A.  RECURSIVE  BAYES  FORMULATION 

The  maximum  amount  of  information  at  time  t.  * in  a pro- 
le 

babilistic  sense*  is  concentrated  in  the  conditional  density 
function  of  the  states  given  all  the  observations  up  to  this 
time*  o ( x^(  k ) ! Zk  ) * where  Zk  represents  the  set  <£(k),  £(k-l)* 
. . . * z_(  0 ) ) . 

By  specifying  a cost  function  one  determines  how  this 
information  is  used  to  obtain  the  best  estimate  [51.  For 
example*  the  minimum  variance  estimate  is  the  mean  of  the 
conditional  density. 

From  the  assumptions  made  in  Section  II*A,  the  random 
processes  x(k)  and  {x(k)*z(k)>  are  Markov  and  the  condi- 
tional density  can  be  written  in  recursive  form  using  Bayes' 
Law*  as  is  shown  in  (51*  {61*  (71,  [81. 


p ( x_(  k ) i Zk  ) = c (k)  .p(x(k)  ! Zk_1)  .p(z  (k)  ! x (k)  ) (3.1) 


where 


p ( x ( k ) !Zk_1) 


p(x(k-l)  ' Zk_1)  .p  ( x ( k ) !x(k-l)).dx(tc-l) 


(3.1a) 


1 / c(k)  = p(z/k)  !Zk-3)  = 

p (jU  k ) : Z1*'3)  .p(z(k)  }_x  (k) ) .dx^k)  (3.1b) 

Given  the  initial  condition  p(x^0)'Z~*)  = p(x(0))  and  the 
probability  law  of  the  system  random  input  and  measurement 
noise#  p(w^k))  and  p(v/k)),one  can  proceed  to  evaluate  these 
equations  and  the  filtering  problem  can  be  regarded  as  hav- 
ing  been  solved  when  p(x^(k)!Z  ) can  be  determined  for  all  k. 

For  any  situation  different  from  the  one  where  the  sys- 
tem and  measurement  equations  are  linear  and  the  apriori 

lc 

distributions  are  Gaussian#  the  density  p(x_(k)!Z  ) is  not 

Gaussian  and  generally  cannot  be  determined  in  closed  form. 

In  the  problem  studied  here#  the  conditional  densities 
in  Eauation  (3.1)  are  nonlinear  functions  of  the  states  and 
the  measurements.  Linearizations#  Taylor  expansions  and 
Gaussian  aoprox i mat i ons  will  be  used#  however#  as  an 
engineering  compromise#  in  order  to  avoid  comolex  numerical 
integrations  in  n dimensions.  This  aDproach  was  chosen  tak- 
ing into  account  the  limitations  on  processing  eauioment 
available  for  the  practical  solution  of  the  addressed  track- 


ing problem 


B.  ONE-STEP  PREDICTION 


Let  us  assume  that  at  time  t^  ^ ? during  a filtering  pro- 
cess? an  approximately  Gaussian  density  p ( x ( k- 1 ) ! Z k-^)  has 
been  obtained  with  mean  value  y Ck-1)  and  covariance  matrix 
M(k-l)  . 

If  one  chooses  the  mean  of  this  conditional  density  as 
the  estimate  at  time  ? 

£(k-l)  = E t x_C  k-1 ) i Zk~'1l  = y (k-1  ) 

the  estimation  error  at  time  k-t  has  the  same  probability 

k— 1 

density  as  o(j^(k-l)|Z  ^ but  with  zero  mean?  i.e.? 

Ete(k-l))  = EKi(k-l)  - x_(  k - 1 ) ) ! 2k~h  = 0 (3.2) 

E (e(k-l)eT (k-1 )J  = P(k-l)  = M(k-l)  (3.3) 

At  time  t,  a new  set  of  measurements  arrives  and  a pred- 
k 

iction  at  that  time  must  first  be  obtained  based  on  the 

measurements  up  to  t,  . 

k-1 

The  dynamics  of  the  system  are  described  by  Eauation 

(2.3)  and  it  can  be  seen  that?  even  if  p ( x ( k- 1 ) ! Zk_1)  were 

i-  1 

really  Gaussian?  the  density  p(x_(k)!Z  ) would  not  normally 

be  because  of  the  non  1 i near i t i es  involved.  Nevertheless  it 

will  be  assumed  that  a reasonably  close  approximation  to  a 

Gaussian  density  exists  and  proceed  to  find  approximate 

values  for  the  first  two  moments  of  the  conditional  one-steo 

. k—  1 

pre-diction  density  p(_x(k)!Z  ). 
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The  first  step  is  to  expand  Equation  (2.3)  in  a Taylor 

series  about  the  estimate  x C k — 1 ) . The  dependence  on  t and  t, 

— k k- 1 

k— 1 

is  dropped  for  simplicity.  Note  that  Z is  known  and  was 
used  in  determining  x_(  k — 1 ) . 

( k ) = f_(  x_(  k — 1 ) ) ♦ g ( x^(  k-  1 ) ) .w(k-l ) = 


= f(5(k-n)  ♦ 

— — 3x 


. Ix(k-l)  - x (k-1 )]  ♦ ...  ♦ 

x(k-l) 


+ a ( x ( k- 1 ) ) . w ( k- 1 ) ♦ ...  = 


jfU(k-l))  - l 
i=l 


3f 

■5x7 


.e  (k-1) 
i 

x(k-l) 


n n 

♦ l l 


- 

3f 

•e.(k-l).e  (k-1) 

3x^3Xj 

i j 

x(k-l) 

i=l  j-1 

♦ q( 5 ( k-1 ) ) • w ( k- 1 ) • 


n 

. I 

i-i 


♦ . . . + 


3g, 

1 

.w  (k-1) 

3x 

i 

x(k-l)  J 

.e(k-l)  ♦ ... 


(3. a) 


where  x^  and  e^  are  individual  components  of  _x  and  (j<  - x_)  # 

g. , is  a column  vector  oi  9 ! and  w^  is  a component  of  w_. 

1 . Linear  Approximation 

In  this  case  only  the  terms  in  the  expansion  which 
are  linear  in  e(k-l)  and  w(k-l)  are  considered#  i.e.# 


x^k)  = _Mx^(k-l))  - <p  (x-1  ) .e(k-l  ) t 
+ g(x(k-l)).w(k-l) 


(3.5) 


where 


f 


i 

S 


' 1 


»(k-l)  = £ Hi<k-l),tk.tk_1)  (3.6) 

From  Equation  (3.5)  the  mean  value  ot  ptjjOcJiZ*  3 ), 
which  is  chosen  as  our  predicted  valuer  is 

£(kSk-l)  = E (_x  ( k ) ! Zk-^I  = f_(_x(k-l))  - 

- $ (k-1 ) .E (e(k-l ) I ♦ g(;(k-l)).Efw{k-l)J 


But  e(k-l)  and  w(k-l)  have  zero  means,  therefore, 

_x  ( k ! k — 1 ) = _f  (_S(k-i)#t  »t  ) (3.7) 

Thus  the  one*step  preaiction  error  has  the  mean 
E te  ( k ! k- 1 ) ] = El(i(k!k-1)  - x ( k ) ) ! £_11  = 

= t[$(k-l).j!(k-l)  - g(_x  (k-1  ) ) .w (k-1  )]  = 0 

The  prediction  error  covariance  matrix,  which  is 
also  the  covariance  matrix  of  p(x^(k)lZ^  ^)  , is,  using  simpli- 
fied notation, 

P(kik-l)  = E re(k!k-l)eT(k!k-l))  = 

= E (k-1)  . 

= E(i)>.e.e_T.  $T-  (j>  . e_.wT.aT  - g . w . e_ T.  cj) T ♦ 9 *T  *°T  ^ 

<16 


e(k-l)  - g(x(k-l)).w(k-l)l  (...) 


T 


From  the  assumptions  about  the  independence  and 


discrete  white  noise  characteristics  of  w*  the  terms 
<p  (k-l).Ete(k-l).wT(k-l)]  .gT(x(k-l))  = 0 
_g(x.(k-l ) ) .E  (w(k-l)  .eT  ( k - 1 )1  . 4>T  ( k- 1 ) = 0 

and 

P(kSk-t)  s <(.  (k-l).P(k-l).  <t.T(k-l)  ♦ 

+ g(£(k-l ) ) .Q(k-i) .gT (x/k-1 ) ) 

2 . First  Order  Approximation 

In  this  development  one  retains  all  terms  in 
tion  (3.4)  up  to  first  order  partial  derivatives,  i ,e 

_x  f W ) s f_(_x(k-l))  “ iji  ( k - 1 ) . e^(  k- 1 ) + 

* q(x(k-l)).w(k-l)  - 
m 

- I A (k-l).w  ( k- 1 ) . e ( k- 1 ) 

, 1 1 “ 

i=l 

where 

A (k-1)  = q.  ( x ( k- 1 ) , t ,t,  ) 

1 3x  — 1 — k k-1 


As  in  the  last  section*  the  predicted  value  is 

_x  ( k i k-1 ) = f(S(k-D)  - * (k-1 ) ,E  le(k-l ))  ♦ 

♦ g(x(k-l)).£Iw(k«l)J  - f A .Efw  (k-l).e(k 

— u , i i — 

i =1 
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(3.8) 


Equa- 


(3.9) 


(3.10) 


-1)3 


From  the  assumptions  about  ^(k-1)  and  e(k-l)*  the 
predi ct i on  is 

_x  ( k } k-1  ) = Hi(k-l)), 

the  same  as  in  Equation  C3.7). 


The  one-step  prediction  error  has  zero  mean  and 
covariance  matrix  given  by*  in  simplified  notation* 

P(kjk-l)  = E fc<(.  (k-1  ) .e(k-l)  - g(£(k-l  ) ) .w(k-l)  + 

+ l Ai(k-l).w  (k-l).eCk-Dl  |s 


r=i 


m 


= E t <().e.eT  . 41 T - <j>  •e»*T*9T  ♦ j>  <J>  .e_.eT  • AT*w  “ 

i-1  1 1 

T T m m ® *T«  T 

• q . w • e • 4>  ♦ g.w*wT#gT  • £ g • w • e1  • A . • w 

i=1  ii 


m 


y A -e.eT.<()T  - y A .e.w1.g1.w  ♦ 

i T i 1 i 


m 


T „T 


i-1 
m m 


i-1 


♦II  *Ai  *wi*wi  1 

i-1  i-1  1 1 J 


-ii 


From  the  assumptions  about  w*  one  has  that 
Ele(k-l).wT(k-l)I  = 0 

E te(k-l ) .eT(k-l  ) .*  (k-1  )1*  = E le  ( k- 1 ) . eT  ( k- 1 ) 1 . 

E [w  (k-1 ) J = 0 
i 

E (w(k-l ) .eT(k-l ) .w  (k-l)l  = EIw(k-l).w  (k-1)]. 

i i 

E leT(  k - 1 ) ] = 0 

«e 


r'w 


i 


. 

*■  - 


C.  UP-UATING  THE  ESTIMATE 

Let  us  assume  that  we  have  an  approximately  Gaussian 
k-1 

density  p(x_(k)  5Z  ) whose  mean  value  is  the  prediction 
x_(k|k-l)  and  whose  second  central  moment  is  the  prediction 
error  covariance  matrix  P(k'k-l). 

A new  set  of  measurements  is  given/  consisting  of  non- 
linear noisy  observations  of  the  states/  represented  by  the 
relationship 

zCk)  = hU(k),tk  ) + v(k) 

It  is  now  necessary  to  process  £(k)  and  extract  the 
information  it  brings.  A new  approximate  Gaussian  density 
p(_x(k)!Z  ) is  assumed  to  be  generated  whose  moments  are  the 
new  estimate  ^(k)  and  the  estimation  error  covariance  matrix 
P(k)  . 

This  new  density  is  given  by  Eguation  (3.1)  or  by 

pU(k)!Zk)  = p(z(k)  U(k)  ) .pU(k)  !Zk-1)  / o(^(k)!Zk_1) 

(3.13) 


k-1 

p (_x  ( k ) ! Z ; is  assumed  to  be  approximately  Gaussian  with 
moments  _a(k!k-l)  and  P(k!k-1). 


Expanding  the  measurement  eauation  around  x(kjk-l)  one 


has 


i(k)  = h(x(k|k-l))  ♦ 


3X 


. Cjx  ( k ) 
x(k| k-1) 


x ( k ! k - 1 ) 1 ♦ ...  t v(k) 
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Retaining  only  the  linear  terms*  one  has 


2.U)  = H_(_x  ( k } k - 1 ) + H ( k ) . lx(k)  - x(k|k-l))  + v(k) 


(3.14) 


where 


H ( k ) = 


h ( x ( k ! k-  1 ) , t.  ) 
— — k 


(3. IS) 


From  F.quation  (3.14)  one  obtains 


E lz(k)  !Zk_1)  = h(x(k|k-l)) 


and  it  is  easily  shown  that 


Var  (z(k)  !Zk_1l  = H(k).P(k|k-l).hF(x)  ♦ R(k) 


Also  from  Eauation  (3.14), 


E C Z_(  k ) ' x_C  k ) 1 S jj(i_(kik-l))  + H(k).(x(k)  - x ( k I k - 1 ) J 


Var  (z ( k ) ! x ( k ) J = R(k) 


If  one  now  assumes  that  p ( ,z  ( k ) | ^-1)  and  p(_z  (k ) ! x ( k ) ) are 
approximately  Gaussian  densities  with  the  moments  qiven 
above,  then  Eauation  (3.13)  becomes 

1/2 

p(x/k)|ZK)  i |H(k ) .P(k : k-1 ) .HT (k  ) ♦ R ( k ) | / 

„ / n 1/2  1/2 

/ [(2ffF/z. |P(kik-l)  | . | R ( k ) | ) . exp  IB) 


where 
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explB)  = exp 


j Ix(k)  - *(kJk-l))T.P"1(k;k-n. 


+ l z ( k ) - h ( x ( k ! k-1 ) ) - H(k).txCk)  - x(k',k-l)]] 


■[« 


"1( k ) . I . . . ] - 


[ Z ( k ) - M x ( k | k“ 1 ))J  • 


[H(k).PCk!k-l).HT(k)  + R ( k ) f1. [ . . . ] 


(3. lb) 


Simplifying  the  notation,  this  exponent  can  be  rear- 
ranged in  the  following  way 

(x  - x)T.r\(x  - x ) + fz  - h ( x ) - H . ( x - x)lT.  R~^.  t . . . ] - 


- [2  - h(i)lT  . [H.P.HT+  RJ_1.  Cz  - h ( x ) 1 = 


= (x  - ; )T  P-^  X - i)  ♦ (x  - ;)THTfrVlCx  - x)  - 


“ 1 2 - h(x)F  R~^H  ( x - x)  - (x  - x)^HTR^[z  - h(x)l  ♦ 


♦ 1 2 - h(i)]T  P_1[2  - h ( X ) } - 


- tz  - h(i)]T  [HPHT  ♦ Rf1lz  - h(i)] 


» T -1  T -1 

= (x  - x)  (P  ♦ H R M) ( x - x)  - 


(x  - i)THTR'1l7  - h ( Q ) 1 - 


lz  - h(x)]TR~Vt(x  - x)  - 


[2  - h(i)JT  [R_1-  (HPHT  t Rf1)  lz  - h ( x' ) J 
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Using  the  matrix  inversion  lemma  19], 


H-1  - (HPHT  +R)"1  = R H ( P~ 1 + HTW1Hr1HTR-1 


Since  R and  P are  symmetric  matrices# 

R-1H(P-1+  HTR-1Hr1Hr  R-1  = AT[P-1  ♦ HrR";1H]A 

where 

A = CP-1  t HTR1H)'1HTR-1 

ana  the  exponent  of  Eauation  (3.16)  becomes 

- T -1  t -1 

(x  - x)  (P  + HR1H)(x  - x)  - 

- (x  - 0THTR-1(Z  - h(i)l  - 

- [z  - h ( x ) 1 R-1H(x  - x)  - 

- [(P_1  + HTfT1HT1HTR"1(2  - h(i))JT[P-1+  HTR"1H][...] 


which  is  a auadratic  form  that  can  be  expressed  as 

(P_1  + hfR^rW1!*  - hc;n]  . (P_1  + HT  R^H]  . 

And  finally  the  conditional  density  of  Equation  (3.16) 
can  be  put  in  the  form 


d ( x ( k ) :zk) 


a ( k ) .exp 


[- 


| lx(k) 


x ( k ) ]T.  P_1(  k ) 


where#  by  definition# 
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K(k)  = i_C  k ! k — l ) ♦ PCk).HTCk).R-l(k).  tz^(k)  - h(x(kik-l),t  )) 

~ — k 

(3.17) 

and 

P_1lk)  = P-1(  k ! k - 1 ) ♦ HT  (k)R_1(k)H(k) 
or , using  the  matrix  inversion  lemma  again, 

P(k)  = P ( k ! k-  1 ) - P(k!k-l)HT(k).(H(k)P(k!k-l)HT(k)  t 

♦ R ( k ) J"1.  H ( k ) P ( k ' k - 1 ) 
or 

P ( k ) = (I  - G(k)H(k)l  .P(kSk-l)  (3.18) 

where 

G(k)  = P(k!k-l).HT(k).[H(k).P(klk-l).hT(k)  + R(k)l 

(3.19) 

From  tquation  (3.19)  one  obtains 
G ( k ) = P(k  ! k-1  ) .HT  (x)  .R_1(k)  . 

(H(k).P(k!k-l).HT(k).R_I(k)  ♦ H"1 

Premult i plying  by  P(k)P  (k),  and  using  the  expression  found 
earlier  for  P~\k)  gives 

G(k)  = P(k).(I  ♦ HTR-1HP  ( k J k-1  ) ] HTR-1. 

[HP(k;k-l)HTR-1  + I)'1  = 
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= P(k)  . CHTR_1  ♦ HTR1HP(k!k-l)HTR_1l  . 

[I  ♦ HP(k  ! k - 1 )HTR-^)-1 


and  final  1 y » 


G ( k ) = PCk)  .HT(k)  .R-^k) 


(3.20) 


Applying  this  result  to  Equation  (3.17)  gives 


i(k)  = i(kik-l)  ♦ G(k).tz(k)  - h(_x(k!k-l),tk)J 


(3.21) 


If  more  terms  were  taken  in  the  expansion  of  _h(_x(k)»tk  ) 
than  the  ones  considered  in  Equation  (3.14)'  different  esti- 
mates and  covariance  matrices  would  be  obtained.  Higher 
moments  of  the  conditional  densities  would  then  be  required. 


we  are  now  in  position  to  restart  the  one-step  predic- 
tion of  Section  III'8  ana  process  a measurement  ocurrinq  at 
time  t , . 


The  initial  steo  in  the  whole  process  is  at  t where  it 

o 

is  assumed  that  an  aoproximate  Gaussian  distribution  exists 
for  the  initial  value  of  the  states'  x^(O).  Choosing  the  ini- 
tial estimate  as  the  mean  value  of  this  distribution'  one 


5(0)  = E (x(0)l 


E le ( 0 ) 1 = E (i(0)  - x ( 0 ) ) = 0 


F 


P(0)  = E[e(0)eA(0)l  = Var[*(0)J 


lhe  set  of  Equations  (3.7 )»  13.8)#  (3.18)»  (3.19)  and 
(3.21)  is  generally  known  as  the  Extended  Kalmam  Filter 
equat ions. 


D.  PARAMETERS  FOR  THE  TRACKING  PROBLEM 


In  order  to  aDPly  the  above  relations  to  the  passive 
tracking  problem#  a series  of  vectors  and  matrices  first 
must  be  determined. 


1.  System  Dynamics 


Applying  definition  (3.6)  to  Equation  (2.10)  yields 


(3.22) 


At  .cos  b ( k ) 


- At.s(k).sin  b(k) 


A t . s i n b ( k ) 


At.s(k).cos  b(k) 


where  At  = ( t,  -t  ) 
k+1  k 


Breaking  up  Equation  (2.11)  into  its  column  vectors 


gives 
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— Ai  .cos  b ( k ) 
2 


g ( x ( k ) * t »t  ) = 
-1  - k+1  k 


1 At 2 . s i n b ( k ) 

2 


At 


(3.23) 


g ( x ( k ) , t # t ) 
— 2 ~ ' k+1  k 


— At2  . s ( k ) . s i n b ( k ) 
2 

•i- At2  . s ( k ) .cos  b(k) 
2 

0 

At 


(3.2a) 


Definition  (3.10)  requires 


A ( k ) = 


- — At2  . s i n b ( k ) 
2 


(3.25) 


0 


0 


0 


—At2  .cos  b ( k ) 


“it 


Ns 


>! 


A 2 ( k ) = (3.26) 


0 

- A t2  .sin  b ( k ) 

- —At2  . s ( k ) .cos 
2 

b(k) 

0 

i At  .COS  b ( k ) 

- - At2  . s(k)  .sin 
2 

b ( k ) 

0 

0 

0 

0 

0 

0 

Assuming  independent  random  forcing  inputs  yields 


Q ( k ) = E (w(k)w  (k)l  = 


a z(  k ) 0 

s 


0 o2U) 
b 


(3.27) 


Equations  (3.8)  and  (3.11)  require  the  product 


g(x(k),t  »t  ) .Q ( k ) ,gT ( x ( k ) , t »t  ) 

- - k+1  k - - k+1  k 


At2. 


o 2(  k ) 


0 

°b<k> 


(3.28) 


where 

a = t At2  . (cos2  b(k)  . 02(k)  ♦ s2  ( k ) . s i n2b  ( k ) .0  2i  < ) ) 
4 s b 

b = -r-  At2  .sin  b ( k ) .cos  b(k).(a2(k)  - s2(k).02(k)) 
4 s b 
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c = ^At  .Isin  b(k).a*(k)  t s4(k).cos  b(k).ob(k)) 
md 

d = Iai.cos  b(k).af(k) 

4 S 

e = - At . s ( k ) .s  i n b(k).ob(k) 

i s ^ At. sin  b(k).a^Ck) 

2 s 

j = yAt.s(k).cos  b(k).a2(k) 


If  the  expanded  state  variable  formulation  is  used# 
the  matrices  qiven  below  are  needed 


Qe(k)  = E lwe(k)*eT(k)l  = 


Q ( k ) | 0 


! 0 


0 0 1 


e n eT 
q .Q  .q 


q.0.gT  | 0 


l o 2 

lAt2 


•Of 


r*  rl 


(3.29) 


(3.30) 


(3.31  ) 


2 . Measurements 


r 


In  the  following  development  Equation  (3.15)  is 

/ 

applied  to  each  of  the  measurements  described  in  Section 
IlfD  to  yield  the  appropriate  linearized  observation  matrix. 


a.  Bearing  Measurement 


From  Equation  (2.2*0  one  has 


_3_  b . = ( x - x )2  / I ( x - x } ♦ ( y - y )2  ] . 

3 l 1 1 1 X 


( ( y - y±  ) / ( x - ) } 


and  then 


■3^  (.£.(  klk-l),^  ) 


(y(kSk-l)  - y^kO/a  l x ( k ! k - 1 ) - x ± ( k ) 1 / aJ  0 0 

(3.32) 


where 


a = d 2(  k ) 


= (x (k ! k-1 ) - 


x^(k))2  + (y(k|k-l)  - y.(k)]2 


(3.33) 


and  di(k)  is  the  predicted  distance  of  the  target  from  buoy 


0.  Frequency  Measurement 


From  Eauation  (2.25)» 
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rfhere  s is  given  in  Equation  (2.26)  and 
i 

£*t  = • C0S(D  * 5 * 

♦ s.sintb  • b ) • b - 
I 

3 

- s.sin(b  - b ^ ) . — b^ 

and  then» 

Hfl(k)  = 4vi(klk*n,v  = 

s UvU  i k-1  )-v  i(k)  1 /£  - ( x Ik  i k-1  1-k^  tk ) 1 /£  6 Y 

where 

f = f (x(kjk-l)»t  ) = f (klk-l).c  / lc  + 
i i - k o 

♦ s ( k ! k-1 ) .cos (b ( k ! k-1 ) - 

b . ( k ) = arc  t an  |jv(k  ! k-  1 ) - y ^ ( k ) 1 / 

(x ( k ! k-  l ) - .t(k),] 

f 2 . s ( k ! k- 1 ) . s i n ( b ( k ! k- 1 ) - 


(3.3a) 


~i  / 
i 

(3.35) 


b ±(  k ) ) ) 

(3.36) 


Y = 


O (k))  / 


(f  (k'.k-D.cl 
o 


f.cos(b(k|k-l)  - b (k))  / (f  (k|k-l).c) 
i i o 


6 ■ v * a in  (3*33) 


c.  Frequency  Difference  Measurement 


From  Equation  (2.2b), 


(f.  - f .)  = f /c  . J-  (S  . - S.  ) 


\ I • I . 

35  1 J 


3?  i J 


where  — s is  qiven  in  Equation  (3.34).  Then, 

35  i 


H (k)  a _L  (f . - f . ) 

di-  3x  l i 


-2-  (f  - f ) = f /c  . (h  h h h ) (3.37) 

3 x i j o 1234 

x(k| k-1) 


where 


a .(y(k!k-l)  - y (k)]  - a .ly(k'k-l)  - y (k)J 

1 i j j 

a ,(x(k!k-l)  - x (k)]  - a .(x(k!k-l)  - x (k)J 

j j i i 

cos(b(k|k-l)  - b (k))  - cos(b(k|k-l)  - b (k)) 

j i 


s ( k ! k- 1 ) . [ s i n ( b ( k ! k- 1 ) - b (k))  - 

i 


sin(b(k|k-l)  - b (k))J 
j 


= s(k!k-l).sin(b(k!k-l)  - 6 (k))  /r 


b (k)  as  in  (3.36)  and  a as  in  (3.33) 
i 
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d.  Frequency  Ratio  Measurement 


From  Equation  (2.2 8)» 


JL  Ct  /f  ) = l/c 
• 1 J 


s ) 
j 


Comparing  with  the  frequency  difference  measure- 
ments one  has 


H (k)  = J-  (f  /f  )|  = 1/f  • H C k ) (3.38) 

ri  3x  i j o di 

x(k| k-1) 


e.  Time  Delay  Measurement 

From  Equation  (2.36)/ 

At  s 1/c  . -L(d  - a.) 

3^  ij  35  i j 

where 

d 2 = ( x - * )2  + (y-y)2 
i i i 


and  then/ 


^ij 


.(  k ) = 


At 

3x  ij 


= l/c  . lh  h 2 0 0) 


(3.39) 


x(kjk-l) 


where 


h L = (x(k',k-t)  - x i ( k ) 1 / dt(k)  - 


PS? 

•s 

< 


li(k!k-l)  - x ^ ( k ) ] / 3.(k) 
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IV.  SPECIAL  CHARACTERISTICS  OF  KALMAN  FILTERS 


i 


* ' 


Equations  (3.7)  and  (3.8)  characterize  the  one-step 
prediction  of  an  Extended  Kalman  Filter;  Equations  (3.18), 
(3.19)  and  (3.21)  represent  the  estimation  update. 

When  implemented  in  a computinq  device,  the  timing  of 
these  steps  is  shown  in  Fig.  5.  The  actual  state  is 
represented  by  a broken  line  and  the  output  of  the  filter  by 
the  solid  line.  In  general,  one  is  dealing  with  vector 
processes  so  Figure  5 is  representative  of  only  one  com- 
ponent. Nevertheless,  this  represent  at i on  is  heloful  in 
visualizing  the  steps  involved  in  the  estimation  process. 

In  a general  problem  the  times  of  occurrence  of  meas- 
urements are  not  equally  separated  ana  are  not  known  in 

advance.  At  a time  t , when  a new  set  of  a measurements 

k-1 

becomes  available,  the  output  of  the  filter  is  still  the 

last  estimate  made  close  to  t , unless  predictions  are  made 

k-2 

at  regular  time  intervals  independent  of  the  spacing  between 
measurement  s . 

A time  period  T is  spent  in  a prediction  phase  to 
P 

evaluate  Equations  (3.7)  and  (3.8),  after  which  the  best 
information  about  the  state  of  the  plant  is  the  predicted 
value  x ( k- 1 ! k-2 ) . 
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Figure  5 : Timing  in  a Kalman  Filtering  Process. 


If  the  times  of  occurrence  of  measurements  were  known 


in  advance/  the  predictions  could  be  made  any  time  the  sys- 
tem was  idle  waiting  for  the  new  measurements  and  T could 


be  made  zero. 


P 


Next/  the  g measurements  are  processed  using  Equations 

(3.18)/  (3.19)  and  (3.21)  and  a time  T is  spent  which  is 

e 

relatively  large  mainly  because  of  a matrix  inversion  of 
dimension  (g  x q)/  in  the  gain  eguat ion/  (3.19).  Only  after 
this  time  has  passed  is  the  new  estimate  x(k-l)  available. 
This  estimate  will  be  the  output  of  the  filter  until  t^  when 
the  process  is  repeateo. 


Some  important  points  in  this  estimation  process  must 
be  stressed: 

---  the  new  estimates  are  only  available  F + F 

P e 

seconds  after  the  arrival  of  new  measurements. 

---  a major  part  of  the  time  period  of  duration  T is 

e 

spent  in  the  gain  equation  due  to  the  required  matrix  inver- 
sion. 

---  from  the  ena  of  the  computation  of  an  estimate 

until  the  arrival  of  a new  set  of  measurements  the  computing 

equipment  is  idle  and  thus  free  to  execute  other  chores. 

---  the  one-step  prediction  is  based  on  a linearization 

of  the  plant  dynamics  over  an  extended  period/  for  example/ 

from  t to  t . 
k-1  k 

---  the  measurement  equations  are  based  on  a lineariza- 
tion about  the  predicted  values. 
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A.  PARTITIONING  OF  MEASUREMENTS 

When  some  of  the  measurements  come  from  statistically 
independent  sources/  steps  can  be  taken  to  alleviate  orob- 
lems  which  arise.  Considering  the  special  case  when  all  the 
measurement  noise  components  are  independent/  the  covariance 
matrix  R(k)  is  diagonal  and  the  following  development 
app lies. 

1 . The  Linear  Case 

Assume  a dynamic  system  with  linear  observations  of 

the  form 

Z_(k)=H(k).>^(k)+v(k)  (4.1) 

where  H(k)  is  not  a function  of  the  states. 

The  estimation  equations  of  Chanter  III  become  the 
standard  Kalman  Filter  equations 

G(k)  = P(k|k”l)HT(k)  (H(lc)P(k!k-l)HT(0  + R ( k ) l-1  (4.2) 

x(k)  = x(k!k-l)  ♦ GOcHz(k)  - H ( k ) .£(  k 1 k- 1 ) 1 (4.3) 

P(k)=lI-G(k)H(k)]P(k!k-l)  (4.4) 

The  measurements  are  assumed  to  occur  simultaneously 
and/  as  shown  for  the  first  time  in  (91/  for  this  linear 
case/  the  results  obtained  by  processing  all  g measurements 
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at  once  using  these  eauations  are  the  same  as  the  final 
results  obtained  by  processing  one  measurement  at  a time. 

In  the  latter  approach,  the  result  of  processing  one 
measurement  component  is  used  in  the  following  computation 
to  process  the  next  measurement  component. 

A general  way  to  show  that  the  two  approaches  yield 
the  same  estimate  after  all  measurement  components  have  been 
processed  is  given  in  the  following  paragraphs. 

In  short  it  is  to  be  proved  that  the  estimates  x(k) 
and  P ( k ) obtained  by  the  use  of  Eauations  (4.2),  (4.3)  and 
(4.4),  when  all  measurements  are  qrouoed  in  a vector  jr,  are 
the  same  as  those  obtained  bv  the  use  of  the  following 
iterative  eauations  a times,  once  for  each  component: 


G . = P.  , . H ,T  / [HP  H?  + R.  1 
i i-I  i i i-1  i i 


(4.5) 


x = X + G ( z - Hx,)  = 
— l -i-I  l—i  i-i-l 


= (I  - G.H  1 x . ♦ G .z. 

l i —i-l  l i 


(4.6) 


(> 

l 

tit 

i 

f 

V *2 

js 

J 


p = 


11  • GiHi]‘Pi-i 


(4.7) 


i - 1,  2,  • • . , a 


where 


= P ( k } k - 1 ) , 


x = x ( k k - 1 ) 
— o — 
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z ( k ) = 


q J 


* ( k ) ♦ 


v 

L q 


and  E (v  1 = 


It  Equation  (4.6)  is  applied  a times#  the  following 


relationship  between  x(k)  (or  x ) and  x(k'k-l)  (or  x ) 

~ - q - -o 


resul t s 


x (k) 


(I  - G H ) [I  - G H 1 ...  T I - G H ] x (k ! k-1 ) ♦ 
q q q-1  q-1  11- 


♦ G z + (I  - GH1G  z + ...  ♦ 
q q q q q-1  q-1 


♦ [I  ” GqHqJtI  ’ • • • < I " G<>  H„  J G,  Z . 


q-I  q-f 


2 2 11 


(4.8) 


and  recursive  use  of  Equation  (4.7)  gives 

P(k)  = (I  - G qH  q 1 (I  - G jH  G1 H1)  P ( k | k-1  ) 

(4.4) 


Comparing  (4.8)  to  (4.3)  and  (4.4)  to  (4.U)  one  sees 
that  the  assertion  can  be  proved  by  showino  the  validity  of 


the  re  1 at  i ons 


II  - G(k)H(k)] 


II  -G  H ] [I 

q q 


G h j 
q-i  q-i 


II  - G H J 
1 1 


(9.  to) 


r (n  x 1)  j 

i 

i 

(n  x 1 ) 

:n  x 1 ) 

G ( k ) = 

i 

II  - GqHq]...[I  - G2H2]Gl  | ... 

1 

1 

(I  - G H JG 

q q q-i 

G 

q 

9.11) 

or*  in  a compressed  way. 


(n  x ( a- 1 ) ) 


( n x 1 ) 


G ( k ) = 


( I - G H 1 G* 
q q 


(a.  12) 


and 


II  - G ( k ) H ( k ) ) = II  - G H ) [I  - G*H*1  (a. 13) 

q q 


* 

where  G is  the  gain  obtained  from  Equation  (3.19)  if  one 
has  a q-1  dimensional  measurement  vector,  and  H* and  R*are 
equivalent  matrices,  as  shown  below 


H ( k ) 


(9. 19) 


R(k) 


(9.15) 
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P(k|k-PH 


(H*P  ( k J It  - 1 ) H 


♦ R] 


(4.16) 


G*  = 


Let  us  expand  Eauation  (4.2)  and  try  to 
it  into  the  form  of  Eauations  (4.10)  - (4.13).  he 


[HPH  + R)  = 


h2phJ 


T 

H PH, 

q 1 


VH2 


W r2 


H.  PHT 

i q 


H2PHq 


H PH1* 

q q 


or 


[HPH  + R)  = 


T 

—3 


—2 


1 


where 


A , = n*PH*T  ♦ R* 


5-2 


H*PHT 


— 3 


T _ 


H PH^  = a* 
q -2 


H PH  + R 

q q q 


manipulate 

have 


V 


(4.17) 


Next  define  the  scalar  c as 


H PHT  ♦ R - H PH  T [H  PH  Tf  R*]_1H  PHT 

q q q q q 

Using  the  definition  of  G*  in  Equation  (4.16)  gives 

H PHT  ♦ R - H G*H*PHT  = 

q q q q q 

: H II  - G*H*]PHT  ♦ R (4.18) 

q q q 


It  is  shown  in  (101  that 

i—  *> 


= A'1  ♦ A-WA-V1 

(4.19a) 

1 1 -2-3  1 

s - A-1  a.  c-1 

(4.19b) 

1—2 

T.-l  -1 
= - a .A.  c 

(4.19c) 

-3  1 

If 

o 

1 

H* 

(4.19a) 

Using  these  re  1 at i onsh i ps , 

• 

PHT  (HPHT  + Rf1  = 


s [PH  B..  ♦ PH  Ab:  ! PH  Ab  ♦ PH1  b ] = 

1 q-3  - 2 q 4 

s [C  ! 01  (4.20) 


Applying  Equations  (4.19a)  - (4.19d)  to  (4.20)/ 

D = - PH*T[H*PH*T  ♦ R*  r1H*PHTc"1+  PH1^1* 

q q 

= - G*H*PHTc_1t  PHTc_]=  (I  - G*H*)PHTc1= 

q q q 

= (I  - G*H*]PHT[H  (I  - G*H*)PHT  + R f1 

q q q q 

Let  P*  = II  - G*  H*  )P  be  the  estimation  error 
covariance  matrix  after  the  processing  of  a-1  measurement 
components/  as  seen  from  Equation  (3.18).  Then/ 

D = P*  HTtH  P*Ht  R T1  = G (4.21) 

q q q q q 

is  the  gain  when  processing  the  qth  measurement. 

Vie  also  have 

C = PH*T[H*PH*T  ♦ R*)-1  + PH*T[H*PH*T  ♦ R *r1H*  PH TH  PH*t. 

q q 

[H*PH*T  + R*  r1c-1  - PHTH  PH*T  [H*PH*t  + R*rLC-1  = 

q q 

s G*  + G*H*PHTH  G*C-1  - PHTH  G*c-1 

q q q q 

Since  c is  a scalar,  one  can  write 

C = G*  - [I  - G*H*)PHTc_1H  G*  = G*  - DH  G* 

q q q 


and.  therefore. 


C s II  • G H ]G 

q q 


(4.22) 


and  f rom  ( 4 . 15 ) # 


■0 


G = |(I  - G H ) G 

q q 


,] 


(4.23) 


GH 


It  is  also  the  case  that 


: (I  - G H ]G*  G 

[_  q q qj 


= (I  - G H )G*H*+  G H 

q q q q 


so 


I - GH  = 


I - G H - (I 

q q 


G H ] G*  H * 

q q 


and#  thus 


I • GH  : [I  • G H J [I  • G*H*I 

q q 


(4.24) 


Equations  (4.23)  and  ( ,m2")  are  the  same  as  Equa- 
tions (4.12)  and  (4.13)#  proving  the  initial  proposition. 

Computationally#  one  can  replace  a (oxq)  matrix 
inversion  and  a (nxq)  x (qxg)  matrix  multiplication  by  only 
q scalar  divisions.  For  an  indication  of  how  much  computinq 
time  can  be  saved#  consider  that  a matrix  multiplication 
requires  approximately  nq2  scalar  mu  I t i o I i c at i ons  anq  addi- 
tions#  and  that  a matrix  inversion  requires  at  least  a mul- 
tiplications and  additions#  not  counting  logic  and  indexing 


75 


- r* 


time  nor  considering  storage  requirements. 

The  partitioning  of  the  measurements  then 

2 

corresponds  to  trading  at  least  q Cn  + q)  multiplications 
and  additions  for  q divisions.  For  a machine  with  average 
multiplication  time  of  10  m i c roseconds , addition  time  of  2 
microseconds  and  division  time  of  12  microseconds,  if  n and 
q are  equal  to  4 a saving  of  about  1.5  ms  is  obtained;  if  n 
and  q are  of  the  order  of  20,  this  saving  could  be  of  tne 
order  of  200  ms  at  each  measurement  time. 

Besides  this  comouting  time  saving,  the  processing 
of  one  measurement  at  a time  provides  partial  estimates  that 
can  be  used  as  improvements  over  the  predicted  values,  as 
shown  in  Figure  6.  Also,  the  final  estimate  is  obtained 
sooner,  although  with  the  same  accuracy  as  when  the  measure- 
ments are  all  processed  simultaneously. 


2.  The  Non  1 inear  Case 


«e  assume  the  same  system  as  before  but  with  non- 
linear observations  of  the  form 

z^  k ) = ( x_(  k ) , t^  ) ♦ v^(k) 


By  linearizing  this  equation  as  shown  in  Chapter 
III,  using  the  definition 

H(k)  = -2-  h ( x ( k | k-1  ) , t ) , 

8x  — — k 

the  proof  of  the  last  section  is  still  valid  for  the 
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k-1) 


\ \+TP 


t +T  +T 
k P e 


measurements 
available . 


predicted  value 
available;  start 
of  processing  all 
q measurements. 


"end  of  processing 
all  q measurements 
estimated  value 
availab le . 


/ \starL  ux  end  of  processing  all  q 

predicted  ^alue  Pr°cessinS  measurements;  final  estimated 
Callable;  start88'0”4  ' value  available, 

of  processing 
first  measurement. 

Figure  6 : Partitioning  of  measurements 
Linear  case. 


linearized  equations 


If#  however#  the  partial  estimates  obtained  by  pro- 
cessing one  measurement  (hopefully  a better  approximation  to 
x_(  k ) than  £(k!k-l)  results)  are  used  to  calculate  the  H 
vector  for  the  processing  of  the  next  measurement#  better 

accuracy  should  be  expected  than  when  all  H vectors  are 

i 

obtained  using  the  relatively  crude  predicted  value.  This 
effect  should  be  particularly  significant  when  the  time 
between  each  group  of  measurements  is  relatively  large. 

Thus#  partitioning  the  measurements  in  an  Extended 
Kalman  Filter  should  provide  improvement  in  estimation  accu- 
racy in  addition  to  the  advantages  pointed  out  in  the  previ- 
ous section.  This  improvement  can  be  visualized  as  shown  in 
Figure  7. 

The  expected  improvement  in  the  estimation  process# 
brought  by  partitioning  the  measurements  in  Extended  Kalman 
Filters#  can  be  explained  in  terms  of  the  conditional  aensi- 
t i es . 


Consider  two  independent  simultaneous  measurements# 


f 

6 


I t 


r 


p(x(k)!Zk)  s pU(k)  !*  (k),z2(k),Zk_])  = 

= p(Zl  (k),z2(k)  !x(k)).p(_x(k)  !Zlt'1)/p(2i  (k),22(k)  !Zk_1)  = 

= lp(i2(k)S*1(k),x(k))  / p(z2  (k)  !Zi  (k),^"1)]  . 

tp(Zl  (k)  !x_(k)).p(_x(k)  :Zk_1)/p(Z;L(k)  |Zk-1)] 

But  t f rom  (3.1) 

p(x(k)  :Zl(k),Zk_L)  = p(Zl  (k)  Jx.(k)  ) .pU(k)  :zk_1)  / 

p ( z ( k ) ! Zk_1) 

also/ 

p ( Z 2 ( k ) ! z^  ( k ) , x_(  k ) ) = p ( z2  ( k ) ! x^(  k ) ) 

because  of  the  indeoendence  assumption. 

And  sor 

o ( x ( k ) i Zk  ) = [ p ( z (k)ix(k))  / p(z  ( k ) ! z ( k ) , Zk_1)  ] . 

- 2 - 2 1 

p(ji(k)|Z  (k)#  Zk_1)  (a. 25) 

The  terms  on  the  right  hand  side  of  (4.25)  are: 

d (j<(  k ) ' z ( k ) , Z11-1)  ---  has  all  the  information  about  the 
states  from  processing  up  to  the  first  measurement  component 
z 1(  k ) . 

o(z2  (k)  Jx_(k))  ---  is  obtained  directly  from  the 
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measurement  equation  and  doesn't  depend  on  the  first  meas- 
urement . 

d ( z 2(  k ) ! z ( k ) # ---  is  obtained  from  the  measurement 
equation  but  using  all  the  results  given  by  the  processing 
of  the  first  measurement  component. 

This  last  term  indicates  that  information  is  lost 
if#  in  a Extended  Kalman  Filter#  the  linearization  for  pro- 
cessing a measurement  is  based  on  the  predicted  value 
instead  of  on  the  value  estimated  after  processing  the 
preceding  measurement  component. 

3 . The  Tracking  Problem 

In  the  tracking  problem  here  addressed#  the  measure- 
ments naturally  occur  at  different  times.  Bearinq  indica- 
tions are  basically  continuous  or#  if  oreorocessed#  succes- 
sive measurements  are  available  within  short  intervals  of 
time.  Frequency  indications  must  be  obtained  from  digital 
processing  of  records  of  considerable  time  length.  Time 
delay  indications  are  also  obtained  from  digital  processing 
but  with  longer  time  intervals. 

In  earlier  work,  [11#  121  anu  (31#  these  measure- 
ments are  forced  to  occur  simultaneously  and  are  processed 
together . 


In  the  orevious  sections  of  this  chapter  it  was 
shown  that#  if  the  measurements  are  all  available  at  tte 
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same  time,  they  should  be  processed  separately  for  better 
computing  efficiency  and,  with  nonlinear  measurements,  for 
possibly  better  estimation  accuracy. 

Since  one  is  now  free  from  having  to  process  the 
measurements  together,  and  encouraged  to  do  so,  why  not  use 
the  idle  time  of  the  computing  equipment  and  process  the 
measurements  at  different  times? 

If  one  does  this  the  decrease  in  computation  time 
obtained  previously  will  be  diminished  because  of  the  extra 
computation  time  required  to  do  many  predictions,  one  for 
each  time  a measurement  is  to  be  processed.  Never t hel ess , 
the  advantages  are  potentially  of  great  value,  that  is, 

---  a measurement  can  be  processed  as  soon  as  it  is 
available,  with  no  waiting  for  other  measurements. 

— measurements  that  occur  more  often  can  be  pro- 
cessed mere  times  than  measurements  ocurring  less  fre- 
quent 1 y . 

---  the  output  of  the  filter  is  updated  frequently. 

— - the  one-step  prediction  is  based  on  lineariza- 
tions through  much  smaller  time-steps  than  before,  hopefully 
with  improved  accuracy  and  fewer  divergence  problems* 

Figure  8 gives  an  idea  of  the  improvement  one 
expects  to  obtain  by  (1)  processing  simultaneous  measure- 
ments separately  and  (2)  processing  the  measurements  *s  they 
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naturally  occur/  when  a Extended  Kalman  Filter  is  being 
used . 


8.  GRAPHICAL  INTERPRETATION 

General  understanding  of  the  actions  taken  by  a Kalman 
Filter  during  the  processing  of  a measurement  can  be 
obtained  through  the  graphical  interpretation  and  visualize* 
t i on  of  these  actions  in  a two-dimensional  problem. 

From  the  conclusions  of  Section  IV#A#  it  is  advantageous 
to  process  any  group  of  measurements  with  statistically 
independent  noise  components  separately.  In  this  section# 
thus#  the  processing  of  only  one  measurement  at  a time  is 
consi dered. 

1 • The  Linear  Measurement  Case 

He  will  consider  a dynamic  system  with  two  state 
variables#  and  a linear  measurement  of  the  form 

z(k)  = H(k).x_(k)  t vCk)  = 

= h^(k).x^Ck)  t h^tkJ.x^tk)  ♦ v(k) 

with  v(k)  a zero-mean#  discrete  white  Gaussian  noise  process 
with  Elv2Ck)J  = r2(k). 

At  time  t^  a prediction  is  made  from  previous  estim 
mates  and  knowledge  of  the  plant  and  system  noise  cherac'* 
teristies*  Equations  (4«3)  and  are  used  for  a linear 


Figure  8 : Comparision  of  processing  policies. 


system,  (5.7)  and  (3.8)  for  a general  nonlinear  system 


The  Kalman  equations  (4.5),  (4.6)  and  (4.7)  are  now 
applied  to  correct  the  predicted  value  according  to  the  new 
measurement . 


Let's  represent  the  quantities  and  relations 

involved  as  shown  in  Figures  9 and  10,  where  e and  e are  a 

“1  “2 

pair  of  orthonormal  basis  vectors  and  a general  vector  is 
represented  by 


Figure  9 : Graphical  interpretation  - Measuremest 
lines  and  gradient. 


M * <»,  * 


g2i1/2 


a = arctan  lh  2^\  * 
6 = arctan  [g  /g  J 


Let's  define  the  function 

<J>  C _><. ) = hixl  * ^2*2  (^.26) 

and  call  "measurement  line"  the  locus  of  points  where  if>  is 
constant.  Some  of  these  lines  are  shown  in  Fig  9, 

Using  this  definition  one  sees  that  the  previously 
defined  Ji  is  the  gradient  of  $ with  respect  to  or  » 

which  in  this  linear  measurement  case  is  not  a function  of 
the  states. 


The  prediction  error  covariance  matrix* 
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a.  The  State  Estimate 

Equation  (4.6)  can  be  represented 
— ~ * 9*lres’dlJa'l 

where 

(residual)  = z(k)  - H ( k ) . i_(  k ! k- 1 ) = 

= 4>  (x_*)  “$(*.')  + v ( k ) 

and  (4.30)  means  that  a vector  correction 
) g j . | res i dua 1 | is  made  to  the  predicted  state 
tion  of  the  vector  n_  which  has  an  angle  g as 
ure  11. 


Equation  (4.5)  gives*  as  shown  in 


P ; 


11  1 


♦ p ' h 

12  2 


p ' h^  ♦ 2p ' h h + p ' h^  + r^ 


11  1 


12  1 2 


22  2 


p ' h ♦ p ' h 

12  1 22  2 


r (4.27) 

(4.28) 

(4.29) 

in  the  form 
(4.30) 

(4.31) 

o f magn i t ude 
in  the  di rec- 
shown  in  F i g- 

Append i * A , 


q 


(4.32) 


(1)  The  Direction  of  Correction 


The  direction  of  g is  seen 

to  be  indepen- 

dent  of 

the 

measurement  noise,  or.  given  h 

and  P(k'k-l)  the 

di rect i on 

of 

correction  is  already  determined  by 

6 = 

arctan 

[(pl2hl  * P22h2J  7 (pUh!  * 

P12h2)J 

(4.33) 

From  Eauations  (4.27)  - (4 

.29)  one  obtains 

p ' 

_ 1 

( a 

- b ) . s i n 20 

(4.34) 

12 

2 

eu 

II 

roll-* 

((a 

+ b)  * (a  - b)  .cos  20  J 

(4.35) 

P22 

_ i 

2 

[(a 

♦ b)  • (a  - b).cos  20] 

(4.36) 

Applying  these  eauations  and  the  definition 

of  a to  (4.33)  gives. 

(a/b  - l).sin2e,  + l(a/b  + 1)  - (a/b  - 1 ) . c os2  0 J . t an  a 

tan  0 = 

(a/b  + 1)  ♦ (a/b  - l).cos20  + (a/b  - 1 ) . s i n2  0.  t an  a 

(4.37) 


which  can  be  reduced  to 

( a/b ).  t ang  . It  an  q . t ana  ♦ 1)  ♦ (tana  • tan0) 
tan  8 = 


( a/b ).( t an  e.  t an  a + 1) 


tan0.[tanoi  - tan0) 


k 


This  equation  shows  that  the  direction  of 
the  correction  generated  by  the  Kalman  Filter#  in  this  two- 
state  problem,  is  a function  only  of  the  direction  ( a ) of 
the  gradient  (fO  and  of  the  alignment  (0 ) and  the  ratio 
(a/b)  of  the  prediction  error  ellipse. 

Three  special  cases  of  interest  occur  when 
(i)  a/b  = 1#  (ii)  b = 0,  and  (iii)  0 = a . 

Case  ( i ) - If  a/b  = 1 the  error  "ellipse" 
is  actually  a circle  meaning  that  the  prediction  has  the 
same  uncertainty  in  any  direction.  No  preferable  direction 
pre-exists  and  the  correction  is  in  the  direction  of  the 
gradient  of  the  measurement  function#  as  can  be  seen  from 
(*1.37) 

tan  0 = (0  + 2. tana)  / (2  ♦ 0)  = tana 
or 

0 = a +•  n tt 


Case  ( i i ) - If  b = 0 the  error  ellipse  is  a 
line  meaning  that  *the  prediction  is  exact  along  the  minor 


axis  and  that  only  corrections  along  the  major  axis  are 
necessary.  As  can  be  anticipated#  the  direction  of  correc- 
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tion  will  be  the  direction  of  the  major  axis#  as  can  be  seen 
by  taking  the  limit  as  b -*•  0 of  Eo.  (4.30)#  which  yields 


tan  0 


t an  0 


or 


6 - 0 ♦ n tt 

For  any  other  intermediate  value  of  a/b  between  1 and  ® , 

the  value  of  6 will  be  between  a and  0 (t  n tt)  . 

Case  ( i i i ) -rthen  0 = a the  gradient  is 

colinear  with  the  major  axis  of  the  prediction  error  ellipse 
and»  f rom  (4.38)* 

tang  = l ( a/b ) . t an  0.  ( t an  20  + 1)  ♦ 0]  / [ ( a /b  ) . ( t an2  0 + 1)  - 0J 
or 

tang  = tan  0 and  0 = 0 ♦ n tt  = a + n u 

Thus  the  direction  of  the  correction  is 
along  the  major  axis  of  the  ellipse  which  is  colinear  with 
the  gradient  of  the  measurement  function. 


(2)  The  Amount  of  Correction 

The  magnitude  of  the  correction 
( | res  i dua  1 ) . 1 g j ) can  be  better  seen  if  one  decomposes  g into 


components  m and  n,  along  the  gradient  and  the  tangent  to 


cos  a = h V j]i  | * sin  a : h / ^ h_ ) 


(4.40 ) 


Using  Equations  (4.32),  (4.39)  and  (4.40), 


p'  h2  + p'  h h ♦ p ' h h ♦ p ' h 2 
11  1 12  1 2 12  1 2 22  2 

Pilhi  * 2d12  hlh2  * 022  h2  * ^ 


1 + Y 


(4.41  ) 


where 


/ (p‘  h2  + 2p  ' h h + o'  h2) 


11  1 


12  1 2 


22  2 


(4.42) 


Consider  some  special  cases.  How  would  the 
new  estimate  be  computed  if  a noise-free  measurement  were 
obtained  and  the  Kalman  fi)ter  were  aware  of  this  fact, 
i . e . , r = 0? 


Note  from  Equations  (4.41)  and  (4.42)  that 
when  the  measurement  is  noise-free,  r = 0,  y = 0 and 


m = 1 / |_hf  , 


a constant. 


No  matter  what  P ( k ! k - 1 ) is,  the  projection 
of  the  correction  into  the  gradient  is  a constant  if  r = 0, 
or,  the  tip  of  the  correction  vector  follows  a line  normal 
to  trie  gradient  when  P(k'k-l)  is  varied. 


Fiqure  12  shows  the  correction  made  to  the 
predicted  state  when  b = 0.  In  Apoendix  C it  is  shown  that 
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for  the  present  case  of  a linear  measurement  function  the 


vector  QA,  in  the  direction  of  h and  length  ( r es i dua 1 1 / | h ^ 
has  its  end  point  A on  the  line 

<(>(*_)  = ( x_*  ) + residual  = 4>  (*. 1 ) + 4>  (_x  * ) + v - $ ( x_'  ) 

or 

( x_)  = ♦ (_**)  ♦ v = z 

which?  for  r = 0 (i.e.  v = 0)  passes  through  the  point  £ = 

x_* . 

Since  the  correction  has  the  direction  of 
the  prediction  error  ellipse  (line)  and  has  OA  as  its  pro- 
jection into  the  vector  h_?  can  be  seen  from  Figure  12 
that  in  this  case  of  b = 0 and  r = 0 the  estimate  will  coin- 
cide with  the  true  state. 

Fiaure  13  shows  the  correction  for  a gen- 
eral error  ellipse  but  still  with  r = 0.  It  is  seen  that  the 
projection  into  h is  the  same  as  before.  The  direction  of 
correction  0 is  somewhere  Detween  a and  0 . 

A general  correction  for  a noisy  measure- 
ment is  shown  in  Figure  14?  with  different  scaling  than  the 
orevious  figures.  The  direction  of  correction  0 is  first 
determined  and  suooose  it  is  as  shown.  The  projection  of  the 
correction  into  _h  would  be  OA  if  r were  zero?  for  r i 0 the 
reduction  factor  is  applied  and  suppose  the  projection  is 
vector  OB  as  shown  in  the  figure.  Takina  a normal  to  _h  from 
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point  d*  the  point  C where  it  cuts  the  direction  of  correc- 
tion determines  the  final  correction  OC  and  the  new  esti- 
mate. 

b.  The  Estimation  Error  Covariance  Matrix 

The  effects  on  the  estimation  error  covariance 
matrix  can  be  seen  in  the  following  way.  From  Equation 
(4.7), 

P C k > = II  - G(k)H(k)]P(k!k-l) 
or 
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Using  the  values  of  g^  and  g ^ from  (4.32)  one 

finds  that 


11= 

th2(pLi  °22 

* -12  > * 

^p2l  ' 

c 

(4.43a) 

• 

Ih  h (o'J 

- P'  P*  ) 

+ p»  r2) 

/ c 

(4.43b) 

12* 

1 2 12 

11  22 

12 

• 

lh2(pr  p' 

- o*2  ) ♦ 

p'  r2l  / 

c 

(4.43c) 

22* 

1 11  22 

12 

22 

where 

c = Ip'  h2  ♦ 2o'  h h ♦ p ' h 2 \ r2  1 (a.  a 3d) 

11  1 12  1 2 22  2 


And  the  estimation  error  ellipse  has  the  orien 


tat  ion 


The  limiting  situations  occur  in  the  special 


cases  when  Ci)  r ■+  ® and  (ii)  r a 0 


Case  ( i ) * When  r is  relatively  large*  little 


information  is  provided  to  the  filter*  the  estimate  is 


essentially  the  same  as  the  prediction*  and  the  estimation 


error  ellipse  is  essentially  the  same  as  the  prediction 


error  ellipse.  This  can  be  seen  from  Equations  (4.43a) 


Case  Cii)  - when  r ~ 0*  all  errors  in  the  pred 


iction,  in  the  direction  of  the  gradient  of  the  measurement 


function*  are  corrected  and  the  estimation  error  ellipse 


becomes  a line  col  inear  with  the  measurement  line.  As  shown 


below*  this  result  is  obtained  from  Equations  (4.43)  and 


(4.44)  when  r is  made  equal  zero 


From  (4.43a)  - (4.43d)»  p .p  is  easily  shown 


when  r = 0*  indicating  that  P(k)  becomes 


to  be  equal  to  P. 


singular.  From  (4.44)* 


= 2(h  /h  ) / Cl  - Ch  /h  )2  J = 

2 1 2 1 

2 

= 2 tan  a / (1  - tan  a)  = tan  2 a 

From  analysis  of  the  individual  terms  of  (<1.430) 
• (4.43d)/  this  result  is  to  be  interpreted  as 

2 6 s 2a  ♦_  tt 
or 

6 s a t_  ir/2 

When  the  minor  axis  of  the  prediction  ellipse  is 
already  zero*  as  in  Figure  12*  i.e.*  o^  = p ' p*  * and  the 
measurement  is  noiseless*  then  P(k)  = 0 and  complete  cer- 
tainty about  the  system  state  is  obtained. 

These  two  last  paragraphs  show  a very  interest- 
ing situation:  if  two  perfect  measurements  are  made  on  a 
two-state  plant  and  the  filter  is  aware  of  this  (i.e.  r = 0 
is  used  in  the  Kalman  equations)*  then  an  exact  estimate  is 
obtained  — the  estimation  error  covariance  matrix  col- 
lapses into  the  zero  matrix. 

For  a genera)  situation  with  normal  measurement 
noise  values*  the  estimation  error  ellipse  is  rotated  toward 
the  measurement  line*  away  from  the  gradient*  with  5 occu- 
pying a position  between  e and  o ♦ v/2. 


r 


... 


pm 


i 


The  general  expressions  for  the  lengths  of  the 
major  and  the  minor  axis  of  the  estimation  error  ellipse  can 
be  obtained  by  the  application  of  Equations  (4.28)  and 
(4.29)  to  (4.43a)  - (4.43c). 

c*  The  Important  Points 

Summarizing  the  principal  points  seen  up  to  now#  j 

---  the  direction  of  correction  in  the  estimate  is 
defined  at  the  end  of  the  prediction  phase  and  doesn't 
depend  on  the  measurements. 

---  if  the  measurement  is  noiseless  and  r = 0 is  used 
in  the  Kalman  filter  equations#  the  projection  of  the 
correction  on  the  gradient  of  the  measurement  function  is 
(residual  / |£|)#  independent  of  P(k{k-1).  The  estimation 

error  covariance  matrix  becomes  singular  indicating#  in  a 
two-dimensional  problem#  that  the  estimation  error  ellipse 
becomes  a line  -~  a point  when  the  prediction  error  ellipse 
is  al ready  a line. 

---  if  the  measurement  is  noisy#  the  projection  of  the 

correction  on  the  gradient  is  multiplied  by  the  reducinq 

2 

factor  (1  / (1  ♦ y The  error  ellipse  is  rotated  toward 

the  measurement  line#  away  from  the  gradient  of  the  measure- 
ment function. 

---  if  the  measurement  is  too  noisy#  the  reducinq  fac- 
tor tends  to  zero  and  the  estimation  is  the  same  as  the 
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prediction. 


2.  The  Nonlinear  Measurement  Case 


Assume  the  same  two-dimensional  system  as  before  but 
with  nonlinear  observations  of  the  form 


z(k)  s h(x(k),t  ) ♦ v(k) 

— K 


a.  The  Extended  Kalman  Filter  Approach 


The  prediction  phase  is  the  same  as  before  but 
the  estimation  equations  are  now  (3.18)*  (3.19)  and  (3.21) 
for  an  Extended  Kalman  Filter. 


Let's  again  define  a measurement  function  by 


$»(  x ) s h(x) 


and  redefine  .h  as  the  gradient  of  0 with  respect  to  x, 


Now  the  components  of  _h  are  no  longer  constants 
but  depend  on  x^.  At  the  predicted  point  these  components 
will  be 


h sih(i(klk-l))  , 


3XX  - 


h =-i_  h(x(k,'k-l)) 

3x2  ~ 


Equations  (9.27)  through  (4.44)  are  all  still 


val id. 


Let's  consider  two  typical  measurement  functions 
and  observe  graphically  how  the  new  estimate  is  obtained. 
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Case  ( i ) 


assume  the  measurement  function  is  of 


the  form 

$(jc)  s a ret  an  Kx^  - a^  ) / ( ^ " a2  71  (4.45) 

The  measurement  lines  will  be  as  shown  in  Figure  15. 

Suppose  now  a noiseless  measurement  is  obtained 
and  is  to  be  processed.  Fig.  16  shows  the  situation  when  the 
prediction  error  ellipse  is  a line  (b  = 0).  In  Appendix  C it 
is  shown  that  for  this  special  measurement  function  the  vec- 
tor OA  of  length  (res i dua 1) / f h f will  end  somewhere  before  a 
point  8 on  the  curve  i|>(jO  s z.  Supposing  its  length  is  as 
shown*  it  can  be  seen  that  an  imperfect  estimate  is 
obtained.  Nevertheless  the  estimation  error  covariance 
matrix*  given  also  by  Equations  (4.431*  will  collapse  into 
the  zero  matrix  indicating*  falsely*  that  a perfect  estimate 
was  obtained. 

Figure  17  shows  a typical  situation  for  a gen- 
eral prediction  error  and  a noise-free  measurement.  The  tip 
of  the  correction  vector  can  be  seen  to  follow  the  line  A 
A'  parallel  to  the  line  <|>(j^)  = <j>(x/)  instead  of  being  along 
the  line  <p(ju)  - <f>(x^).  As  seen  in  the  previous  section* 
the  estimation  error  ellipse  will  erroneously  be  determined 
as  a line  also  along  A - A'*  indicating  falsely  the  direc- 
tion of  the  true  state. 

Case  ( i i ) - assume  the  measurement  function  is 


of  the  form 


residual 


Figure  17  : Nonlinear  measurement  I - General 
correction  for  r=0. 


♦ (*„ 


(4.46) 


The  measurement  lines  are  as  shown  in  Figure  18,  This  same 
figure  shows  the  new  estimate  obtained  when  the  measurement 
is  noiseless  (r  = 0)  and  the  prediction  error  ellipse  is  a 
line  (b  = 0).  The  same  observations  of  Case  (i)  are  valid. 
As  shown  in  Appendix  C the  point  A is  on  the  curve  <$,(x)  = z 
for  this  special  measurement  function. 


These  figures  have  their  values  mostly  because 
they  show*  in  a simple  way#  the  large  errors  that  can  be 
made  in  the  processing  of  nonlinear  measurements  by  an 
Extended  Kalman  Filter.  Two  types  of  errors  are  generated: 
one  in  the  determination  of  the  new  state  estimate#  the 
other  is  seen  in  the  fact  that  the  estimation  error  covari- 
ance matrix  P becomes  a poor  approximation  to  the  true  error 
covariance#  making  the  filter  non -optimal.  These  difficul- 
ties are  a function  of  the  nonlinearities  involved  ana# 
apparently#  are  more  significant  when  the  predicted  values 
are  poor  and  when  the  measurement  noise  is  assumed  to  be 
small  when  evaluating  the  extended  Kalman  filter  equations. 


b.  The  Iterative  Approach 

A way  of  correcting#  in  part#  these  deficiencies 
is  the  adoption  of  iterative  procedures.  The  first  improve- 
ments can  be  obtained  if#  after  the  Extended  Kalman  Filter 
has  been  applied  to  correct  the  prediction#  the  filtering 
process  is  repeated  but  with  the  1 i near i zat i ons  made  about 
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the  just  obtained  estimate*  hopefully  a better  approximation 


to  the  true  state  than  the  crude  prediction. 

Analytically  the  process  consists  of  replacing 
the  estimation  equations  (3.18)*  (3.19)  and  (3.21)  by  the 
iterative  equations 

/ 

x = x ( k S k-1 ) ♦ G Iz(k)  - h(x(k'k-l)>l  (4.47) 

i+1  1 

G1  = P(k!k-l)Hj(H1  R(kik-l)^  ♦ R)"1  (a. 48) 


_*1+1  = *'  ♦ g^Cz  - *(5^)  ♦ - <frO<')] 

or 

* x.'  ♦ ai«  (z  - <|>  (_x  ± J J - gi.t^(x_' ) - ip  (x^)J  (4.52) 

These  equations  can  be  applied  until  there  is  no 
significant  difference  between  consecutive  estimates#  and  2 
or  3 iterations  are  normally  enough. 

The  effect  of  these  equations  can  be  easily 

visualized  for  the  simple  case  of  b = 0 and  r = 0.  The 

sequence  of  Figures  14  « 21  show  this  effect  for  one  of  the 

special  measurement  functions  already  studied.  In  Figure  14 

the  first  estimate  is  obtained  as  done  previously  in  Figure 

lb.  The  vector  0A,  which  is  smaller  than  OB  as  shown  in 

Appendix  C»  represents  the  projection  of  the  correction  on 

the  gradient  h . A normal  to  h from  point  A cuts  the  ored- 
—q  — o 

iction  error  line  (b  = 0)  at  the  first  estimate  • 

In  Figure  20  the  terms  of  Equation  (4.52)  are 
shown  for  i - 1.  Two  correction  terms  are  applied  to  the 
prediction  x_* , both  in  the  direction  of  g_  but  with  dif- 
ferent residuals.  The  direction  of  ’s  again  along  the 
prediction  error  line  (b  - 0). 

The  gradient  of  the  measurement  function  at  the 
first  estimate  is  _h^  . The  projection  of  the  first  correc- 
tion (g^Iz  “ ( *_  i)  J ) into  h_j,  is  smaller  than  CE#  say  CO. 

Taking  the  normal  at  D one  gets  the  first  correction  term 


1 1 1 


CH 


The  projection  of  the  second  correction  term* 
(-c^  [ <J>  (x_* ) • + into  Jh  ^ is  smal  ler  than  CG*  say  CF. 

The  normal  from  point  F determines  the  second  correction 
term  (IC). 

Adding  the  three  vectors  CH  and  IC  one  has 

the  second  estimate*  shown  in  Figure  21*  and  the  next  itera- 
tion can  be  processed. 

If  in  Equation  (4.47)  one  modifies  the  residual 
by 

z ( k ) - h ( x ( k ! k-1 ) ) = z(k)  - lh(x  ) ♦ 

“ ~ ~ i 

♦ H^.(£(k'k-1)  - x^)  + ...]  = 

s'  z(k)  * h ( i ) - H.(x(k}k-l)  - xl 
- — — i.  i — — i. 

then  Equation  (4.47)  becomes 

x = i(k!k-l)  + G (z ( k ) - h(i  ) - H lx (k i k-1 ) - 5 )) 

- i+1  — i ~ i i“  ~ 1 

(4.53) 


This  variation  in  t'he  residual  is  very  small  to 
affect  considerably  the  estimates  but  in  till  it  is  shown 
that  its  use  produces  a better  agreement  between  the  calcu- 
lated estimation  error  covariance  matrix  and  the  true  error 
covariance*  when  processing  a noisy  measurement. 
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3.  Summary 


The  direct  application  of  these  graphical  interpre- 
tations to  a general  problem  is  impractical.  Nevertheless 
the  basic  points  observed  here  can  be  generalized  for  Kalman 
Fi 1 ters: 


---  the  "direction"  of  correction  in  the  estimation 
is  defined  by  the  predicted  quantities  and  doesn't  depend  on 
the  measurements. 


---  the  amount  of  correction  in  an  estimation,  along 
the  direction  defined*  is  the  result  of  a weighting  process 
between  the  uncertainty  in  the  measurement  and  the  confi- 
dence in  the  predicted  values.  This  correction  is  a function 

2 

of  the  factor  (1  / (1  )]  where*  in  general*  for  a n- 

dimensional  system 


f2=  ,!  'ii  j-i  °ijVj 


(4.51) 


---  the  errors  in  the  estimate  from  the  processing 
of  a nonlinear  measurement  can  be  large  and*  most  important* 
the  estimation  error  covariance  matrix  may  become  a-  very 
poor  approximation  of  the  true  error  covariance*  forcing  the 
filter  to  inaccurately  weight  future  measurements.  This 
problem  depends  on  the  nonlinearities  involved  and  seems  to 
be  worse  when  the  predicted  values  are  poor  and/or  the  meas- 
urements are  considered  by  the  filter  to  have  relatively  low 


noise. 
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---  the  errors  indicated  above  may  be  reduced  with 
the  use  of  the  iterative  equations  (4,48)#  (4.49)  and  (4.47) 
or  (4.53)*  instead  of  the  standard  approach  given  by  the 
Equations  ( 3. 18) * ( 3 . 19)  and  (3.21)*  at  the  expense  of 
greater  computing  effort. 
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V.  COMPUTER  SIMULATION 


A.  GENERAL  CONSIDERATIONS 

A computer  simulation  was  implemented  in  a Unix  Operat- 
ing System  running  on  a POP  11/50  with  the  FP  11B  floating 
point  processor,  located  at  room  506  of  Soanagel  Hall,  NPS. 
The  whole  program  is  interactive  and  is  expected  to  be 
self-explanatory  to  a user.  Its  block  diagrams  and  coding 
are  presented  in  Appendix  E. 

The  graphical  outputs  were  generated  on  a Tektronix 
401A-1  terminal.  The  presentation  of  the  results  of  the 
various  simulations  is  simplified  by  the  adoption  of  the 
following  conventions. 

The  metric  system  was  used  in  all  calculations  and  out- 
puts. The  axes  are  dimensioned  in  kilometers.  Heading  and 
bearing  angles,  expressed  in  radians,  are  measured  counter- 
clockwise from  the  positive  X-axis. 

Figure  22  shows  a typical  plot.  Points  of  the  true  track 
are  represented  by  the  sysmbol  "t"  and  interconnected  by  a 
solid  line.  Buoys  are  represented  by  a distinctive  boxed 
symbol.  The  portion  of  the  true  track  that  is  in  the  range 
of  a buoy  is  bounded,  if  necessary,  by  dash-dotted  lines. 


The  signals  received  by  the  sonobuoys  are  assumed  to  be 
sent  to  a ship  or  aircraft  where  they  are  processed  by  spe- 
cialized equipment  to  generate  certain  measurements.  The 
major  characteristics  of  these  measurements  are  printed  on 
top  of  each  figure  and  obey  a code  whose  general  form  is 

L N 1 N2  N5  S 

where 

L - a capital  letter  indicating  a type  of  measurement. 
B stands  for  bearing;  F for  frequency?  D for  frequency 
difference?  R for  frequency  ratio?  and  T for  time  delay. 

N 1 - an  optional  number  indicating?  when  necessary?  the 
buoy  that  generates  the  measurement.  If  this  number  has  two 
digits  it  identifies  two  buoys?  thus?  T 12  indicates  a time 
delay  measurement  between  buoys  1 and  2. 

N2  - a number  indicating  how  many  seconds  after  the 
start  of  the  problem  this  measurement  will  be  first 
obtained.  It  is  also  optional. 

N3  - a number  indicating  the  period  between  consecutive 
measurements . 

S - a number  indicating  the  standard  deviation  of  the 
simulated  zero-mean  Gaussian  noise  that  is  associated  with 
this  measurement?  followed  by  an  ootional  small  letter  indi- 
cating dimension:  d for  degrees?  h for  Hertz  and  s for 
seconds . 
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The  interpretation  of  the  code  presented  in  Figure  22  is 
then  that  both  buoys  l and  2 provide  bearing  and  frequency 
measurements  with  a period  of  100  seconds  between  consecu- 
tive set  of  measurements.  Buoy  1 provides  its  first  set  of 
measurements  100  seconds  after  the  start  of  the  problem; 
buoy  2*  perhaps  because  of  its  limited  range»  provides  its 
first  set  of  measurements  only  after  500  seconds  of  the 
start  of  the  problem.  The  standard  deviation  of  the  measure- 
ment noise  is  5 degrees  for  the  bearings  and  0.04  Hertz  for 
the  frequencies. 

The  filter  receives  the  measurements  for  processing.  It 
is  given  some  initial  conditions  and  the  initial  position  is 
marked  in  the  figures  by  the  first  diamond  with  the  letters 
I.P.  on  top. 

As  the  target  moves  ana  the  measurements  are  simulated* 
the  filter  generates  estimates  of  the  target  parameters 
which  are  represented  by  small  diamonds*  interconnected  by 
dashed  lines  to  represent  the  estimated  path.  The  diamonds 
and  the  symbols  correspond  to  the  same  time  points. 

Associated  with  each  estimated  point  the  filter  provides 
an  estimation  error  covariance  matrix  which  indicates  the 
amount  of  confidence  or*  conversely*  uncertainty  that  should 
be  given  to  this  estimate.  An  interesting  way  to  show  this 
in  a trajectory  plot  like  Figure  22  is  by  plotting  the  one- 
sigma  error  ellipses  associated  with  each  estimated  point. 
These  error  ellipses  have  their  orientation  and  axis  dimen- 
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sions  given  by  Equations  (4,27)  - (4,29),  and  can  be  seen  in 
Figure  23. 

These  figures  were  obtained  by  running  the  simulation 
once.  They  may  represent  a typical  real-time  result. 
Nevertheless  they  do  not  give  the  expected  or  average 
behavior  of  the  filter  to  the  situation  simulated.  To  obtain 
this  average  behavior  many  Monte  Carlo  runs  have  to  be  exe- 
cuted and  their  individual  results  processed  to  generate  the 
sample  statistics.  The  sample  statistics  formulas  for 
obtaining  the  sample  means  and  covariances  can  be  found,  for 
example,  in  page  86  of  Ref.  (21. 

In  this  work  all  the  statistical  results  were  obtained 
from  200  MonteCarlo  runs.  Figure  24  shows  the  average 
behavior  of  the  filter  after  these  runs.  Now  the  small  dia- 
monds represent  the  sample  mean  estimated  position  and  the 
ellipses  show  the  spreading  of  the  results  about  the  mean 
values.  Note  that  at  the  initial  position  point,  since  in 
all  runs  it  was  the  same,  the  ellipse  is  the  point  itself. 

To  observe  the  errors  in  the  estimation  of  the  speed  and 
heading  of  the  target  different  plots  have  to  be  used.  Fig- 
ure 25  shows  the  sample  mean  error  in  the  estimation  of 
speed  as  a function  of  time.  Figure  26  shows  the  sample 
mean  error  in  the  estimation  of  the  heading  of  the  target. 

Other  types  of  plots  or  outputs  could  be  generated  from 
the  results  given  by  the  program  simulation,  like  the  sample 
covariance  between  estimated  velocity  and  bearing,  if  it 


122 


»f>  U>  Ck> 


a.ee  1.50  3.00  4.50  6.00  7.50  9.00  10.50  12.00 

(sec)  (*100 


Fi3ure  26  : Error  .in  .heading  .estimation . - 
Sample  means. 


ever  becomes  necessary 


B.  SELECTED  RESULTS 

As  has  been  seen  in  this  work,  the  number  of  variables 
in  the  passive  tracking  problem  is  very  large.  Some  specify 
the  true  target  track  and  the  true  position  and  characterise 
tics  of  the  sonobuoys;  others  define  the  measurement 
scheduling  and  noise  characteristics.  The  nonlinear  filter 
that  processes  the  measurements  has  its  own  variables  and 
must  also  assume  values  for  the  other  variables,  which  may 
be  different  from  the  true  ones. 

For  any  combination  of  these  variables  a new  problem  is 
formed  and  its  simulation  may  generate  a number  of  interest- 
ing conclusions.  The  results  presented  are  representative  of 
those  obtained  f rom  the  execution  of  the  simulation  program 
described  in  Appendix  E.  In  most  of  the  results  presented  a 
certain  reasonable  initial  condition  for  the  filter  was 
assumed;  the  position  of  the  buoys  was  assumed  constant  and 
known;  the  ranges  within  which  reliable  bearing,  frequency 
and  time  delay  measurements  could  be  obtained  yere  assumed 
equal  . 

1 . Non-maneuvering  Target 

a.  One  Buoy  in  Action  at  a Time 

Figure  2 7 shows  the  behavior  of  the  filter  in  a 
situation  where  only  one  buoy  is  in  contact  with  the  target 
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at  any  instant.  The  range  and  position  of  the  buoys  were 
pre-arranged  to  satisfy  this  reauirement.  Target  speed  is 

5.0  m/si  its  headinq  is  315  degrees  measured  from  the  posi- 
tive X-axis.  The  initial  assumptions  of  the  filter  include 

6.0  m/s  for  speed  and  325  degrees  for  heading.  The  measure- 
ments were  processed  as  suggested  in  121.  Frequency  and 
bearing  measurements  are  available  every  100  seconds  and 
processed  simultaneously;  standard  deviations  of  5 degrees 
and  0.04  Hertz  were  assumed  for  the  noise  components. 

The  frequency  measurement  provides*  indirectly* 
some  range  information  and*  this  way*  complements  the  bear- 
ing measurement.  This  range  information  can  also  be  obtained 
by  judiciously  positioning  the  next  buoys  or  by  having  two 
or  more  buovs  in  contact  with  the  target  at  one  time. 

In  Figure  28  the  frequency  measurement  was  elim- 
inated and  the  simple  bearing  measurement  allowed  the  esti- 
mator to  track  the  non-maneuver i ng  target  almost  as  well  as 
before*  because  of  the  special  position  of  the  buoys.  Note 
that  the  ellipses  tend  to  align  with  the  measurement  lines. 
This  can  be  explained  from  the  discussions  of  Section  IV*B 
and  Appendix  D. 

There  are  many  ways  to  improve  this  tracking.  A 
very  small  improvement  can  be  obtained  if  less  noisy  or  more 
frequent  frequency  measurements  could  be  obtained.  To  get 
more  precise  frequency  indication  one  needs  better  resolu- 
tion from  the  F FT  processors;  to  obtain  better  resolution 
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27  : .True  and  estimated  trajectories  for  a 
non-maneuvering  target  - Simultaneous 
frequency  and  bearing  measurements. 
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one  needs  longer  records  which*  in  turn*  introduce  more 
errors  due  to  t arget -sensor  geometry  variations  durinq 
record  time.  To  reduce  the  period  between  consecutive 
independent  frequency  indications  one  has  to  sacrifice  reso- 
lution but  this  would  reduce  even  more  the  effect  of  the 
frequency  measurement  on  tracking  accuracy.  Besides*  an 
emitted  frequency  of  300  Hertz  from  a target  with  a relative 
velocity  of  5 knots  toward  a buoy  is  shifted  by  only  about 
0.5  Hertz  and  it  can  be  seen  that  not  much  resolution  can  be 
spared. 

The  bearing  information  is  basically  continuous 
and  it  seems  reasriable  to  admit  that  only  economical  limi- 
tations exist  for  obtaining  less  noisy  or,  equivalently* 
more  frequent  bearing  measurements.  Figure  29  shows  the 
effect  of  processing  only  bearing  measurement s * but  more 
frequently  (at  25-second  intervals).  In  Figure  30  the  stan- 
dard deviation  of  the  measurement  noise  was  reduced  to  1 
degree.  Note  that  the  alignment  of  the  ellipses  with  the 
measurement  lines  is  more  pronounced  with  less  noisy  meas- 
urements. Note  also  that  the  error  in  range  while  the  first 
buoy  is  i-n  contact  with  the  target  is  not  corrected  by 
improving  the  accuracy  of  the  bearing  measurements. 

Figures  31  and  32  show  the  mean  errors  in 
estimating  the  velocity  and  bearing  of  the  target.  The  sym- 
bol is  associated  with  Fig  27,  the  triangles  with  Fig  28 
and  the  diamonds  with  Fig  29. 
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Figure  29  : True  and  estimated, trajectories  for  a_, 
non-maneuvering  target  - Reduced  inter 
val  between  consecutive  measurements. 


Since  the  contribution  of  the  frequency  measure- 
ments is  so  small  compared  to  that  of  the  bearing  indica- 
tions in  this  case*  the  partitioning  of  the  measurements  has 
very  little  effect  on  tracking  accuracy.  Also  the  iterative 
techniques  suggested  in  Section  IV/B/2  have  no  marked  effect 
on  the  processing  of  the  frequency  measurements  for  the  same 
reason . 


b.  Two  Buoys  in  Action  at  a Time 

The  placement  of  another  buoy  in  action  can 
greatly  improve  the  tracking  accuracy/  as  shown  in  the  next 
figures.  The  buoys  are  now  arranged  so  that  two  buoys* 
instead  of  one/  are  in  contact  with  the  target  at  any  time. 
The  scaling  is  changed  so  that  the  estimation  errors  may  be 
better  appreciated. 

In  Figure  33  the  added  buoy  is  a Lofar  buoy  that 
can  only  provide  frequency  measurements.  The  two  frequency 
and  one  bearing  measurements  available  every  100  seconds 
were  processed  simultaneously.  Note  the  improvement  in  accu- 
racy when  buoys  3 and  d gain  contact  with  the  target  in 
replacement  of  buoys  1 and  2.  This  also  is  explained  from 
the  discussions  of  Section  IV/6. 

In  Fiqure  3d  all  buoys  contribute  with  frequency 
and  bearing  measurements  which  are  still  processed  simul- 
taneously every  100  seconds.  The  position  of  the  buoys/  in  a 
line  perpendicular  to  the  true  track/  is  not  a very  good 
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Figure  3,3 . : . True  and.  estimated  trajectories  for  a 
non-maneuvering  target  - Simultaneous 
processing  of  one  bearing  and  two  fre 
quency  measurements. 


position  for  two  Difar  buoys 


Nevertheless  the  resulting 


tracking  accuracy  is  good  for  this  non-maneuver i ng  target. 

In  Figure  35  the  two  frequency  measurements/ 
which  require  a five-state  plant  for  their  processing/  were 
replaced  by  a frequency  difference  measurement  with  0.04 
Hertz  of  standard  deviation  in  its  noise  component/  which 
requires  only  a four-state  plant/  as  discussed  in  Section 
II/D/2.  The  results  are  basically  the  same/  mainly  because 
they  are  mostly  dependent  on  the  bearing  measurement  and  not 
on  the  frequency  information.  The  use  of  the  frequency  ratio 
measurement  also  gave  the  same  kind  of  results  when  a very 
low  noise  measurement/  with  about  0.0001  units  of  standard 
deviation/  was  considered. 

In  Figure  36  the  frequency  measurements  were 
eliminated  and  only  bearing  measurements  were  processed. 

Because  of  the  bad  position  of  the  buoys/  on  a line  oerpen- 

' 

dicular  to  the  true  track/  the  accuracy  of  the  tracking 
along  that  line  is  not  good.  The  ellipses  show  a reasonable 
spreading  of  the  estimates  along  perpendi cul ars  to  the  true 
track.  Figure  37  shows  the  same  plot/  with  different  scal- 
ing/ for  the  first  1000  seconds  of  the  oath. 

For  different  positions  of  the  buoys/  especially 
if  one  moves  the  buoys  off  the  perpendi cular  to  the  true 
track/  better  tracking  can  be  obtained  as  shown  in  Figure 
38.  The  tendency  of  the  ellipses  to  align  with  the  measure- 
ment lines  explains  the  improvement  obtained. 
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By  reducing  the  period  between  consecutive  bear- 
ing measurements  and  processing  them  separately  one  may 
obtain  the  accuracy  shown  in  Figure  39.  The  ellipses  are 
like  the  ones  of  Figure  38  but  with  reduced  dimensions. 

Good  results  can  also  be  obtained  if  one 
i nc  1 udes  the  time  delay  measurements . Figure  <10  shows  the 
result  when  they  are  applied  to  this  situation.  Notice  that 
the  ellipses  are  very  small  reflecting  the  fact  that  very 
good  information  about  the  position  of  the  target  is  aiven 
by  a low-noise  time  delay  measurement.  Also  the  ellipses 
tend  to  align  with  the  tangent  to  the  measurement  lines 
defined  in  Section  IV, B which,  in  this  case  of  time  delay 
measurements,  are  hyperbolas. 

2 . Maneuvering  Target 

The  next  paragraphs  discuss  results  obtained  when  a 
selected  target  maneuver  is  simulated  ---  a total  turn  of 
180  degrees  at  9 degrees  per  minute  with  a constant  speed  of 
5.0  m/s . 


a.  One  Buoy  in  Action  at  a Time 

Figure  41  shows  the  tracking  obtained  by  pro- 
cessing freguency  and  bearing  measurements  provided  by  a 
single  buoy  located  close  to  the  center  of  curvature  of  the 
target  path.  The  filter  takes  some  time  to  react  to  the 
target  maneuver  and,  when  it  does,  an  over-correction  is 
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applied.  The  result  is  that  the  estimated  path  is  very  dif- 
ferent from  the  true  path  of  the  target.  About  the  same 
result  is  obtained  when  the  buoy  is  placed  on  the  other  side 
of  the  track,  i.e.,  a small  delay  to  react  to  the  maneuver 
followed  by  an  over-correct i on. 

In  Figure  42  the  frequency  measurements  were 
eliminated  and  only  bearing  indications  were  processed  to 
generate  an  estimated  track  very  close  to  the  one  of  Figure 
41.  The  alignment  of  the  ellipses  with  the  measurement 
lines  is  again  very  clear  and  may  suggest  where  two  buoys 
should  be  placed  in  order  to  improve  the  tracking. 

b.  Two  Buoys  in  Action  at  a Time 

Figure  43  shows  what  can  be  obtained  by  process- 
ing the  measurements  provided  by  two  buoys.  The  range  of  the 
buoys  was  adjusted  so  that  both  are  in  contact  with  the  tar- 
get during  all  the  path  shown.  Both  are  Difar  buoys  and  pro- 
vide bearing  and  frequency  measurements  with  a period  of  100 
seconds,  which  are  processed  simultaneously  by  the  filter. 

In  Figure  44  the  measurements  were  processed 
separately  and  a small  improvement  in  the  tracking  was 
obtained  at  the  end  of  the  path,  although  during  the 
maneuver  the  simultaneous  processing  worked  better.  It  was 
observed  in  other  simulations  that  the  degree  of  improvement 
obtained  by  partitioning  the  bearing  and  the  time  delay 
measurements  is  dependent  on  the  maneuver  of  the  target  and 


the  relative  position  of  the  buoys.  The  worst  situation  is 


when  the  target  maneuvers  toward  the  buoys.  This  happens 
when  the  partial  estimates  are  worse  than  the  predicted 
values  and  occur  sometimes  during  fast  maneuvers*  or  slow 
maneuvers  with  long  intervals  between  consecutive  measure- 
ments . 

In  Figure  **5  the  low-noise  time  delay  measure- 
ment was  added  and*  instead  of  reducing  the  errors*  the 
tracking  became  worse.  The  direction  of  the  ellipses  is 
clearly  along  the  tangent  to  the  hyperbolas  of  position, 
indicating  that  the  filter  is  uncertain  about  the  range  of 
the  target  to  the  buoys  during  all  the  maneuver.  This 
characteristic  of  the  measurement  lines  in  this  case*  and 
the  fact  that  the  target  is  modelled  by  the  filter  as  fol- 
lowing basically  a straight  path,  explain  the  form  of  the 
estimated  track.  As  soon  as  the  target  reassumes  a constant 
path  the  filter  starts  to  recover. 

In  Figures  46  and  47  the  iterative  techniques 
described  in  Section  IV*B,2*b  were  applied  to  the  frequency 
and  time  delay  measurement s . In  Figure  46  only  one  itera- 
tion was  executed?  in  Figure  47  the  iterated  formulas  were 
applied  three  times.  The  improvement  is  clearly  seen  by  com- 
paring with  Figure  45. 

As  with  the  partitioning  of  the  measurements,  it 
was  observed  from  other  similar  simulations  that  the  appli- 
cation of  the  iterative  techniques  does  not  always  provide 
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Figure  45  : True  and  estimated  trajectories,  for  a 
maneuvering  target  - Separate  proce- 
ssing of  bearing,  frequency  and  time 
delay  measurements. 
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..Figure.  46,.:.  True,  and  estimated. . traj.ectpr.ies.  fo.r  a 
maneuvering  target  - One  iteration  on 
frequency  and  time  delay  measurements. 
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Figure  47  r_Irue  .and.  estimated  .trajectories  f or.  a .... 

maneuvering  target  - Three  iterations 
on  frequency  and  time  delay  measurements 


better  tracking.  It  depends  on  the  type  and  direction  of  the 
maneuvers*  the  position  of  the  buoys  and*  consequently*  on 
the  position  of  the  measurement  lines. 

In  the  case  of  the  non-maneuver i ng  target  the 
use  of  two  buoys  instead  of  only  one  was  almost  always  suf- 
ficient to  provide  good  tracking.  If  for  the  maneuvering 
target  one  uses  three  buoys  judiciously  positioned*  instead 
of  only  one  or  two*  one  can  also  obtain  good  tracking  as 
clearly  shown  in  Figure  48. 

C.  SUMMARY 

Below  are  listed  the  most  important  facts  observed  from 
the  simulations: 

---  The  tendency  of  the  error  ellipses  to  align  with  the 
measurement  lines  was  observed  in  Section  IV*B.  The  simula- 
tions show  that  the  sample  covariance  ellipses  also  follow 
this  trend. 

---  The  tendency  of  the  error  ellipses  to  align  with  the 
measurement  lines  can  be  of  great  help  in  the  practical 
solution  of  filtering  situations.  If*  for  example*  the  error 
ellipse  is  very  thin*  i.e.*  have  a very  high  ratio  a/b*  then 
the  best  measurement  to  process  is  the  one  whose  measurement 
line  is  approximately  perpendicular  to  the  major  axis  of  the 
ellipse.  In  the  tracking  problem  this  can  be  obtained  by  a 
bearing  measurement  from  a buoy  placed  along  the  perpendicu- 
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Figure  48. True  and  estimated,  trajectories  for  a. 

maneuvering  target  - The  influence  of 
a third  buoy. 


IX  IX 


Jar  to  the  major  axis  of  the  ellipse*  or  by  a time  delay 
measurement  from  two  buoys  placed  along  the  direction  of  the 
major  axis  of  the  ellipse. 

---  The  frequency  measurement  have  a very  small  effect 
on  the  tracking  accuracy. 

---  The  use  of  frequency  difference  or  frequency  ratio 
measurements  improve  computing  efficiency  (by  eliminating 
the  need  for  having  the  rest  frequency  as  a state  in  the 
filter  model)  without  noticeable  effect  on  accuracy. 

---  The  partitioning  of  the  measurements  improves  com- 
puting efficiency*  may  provide  better  tracking  accuracy  and 
allows  the  processing  of  measurements  as  they  occur,  and 
thus  is  of  great  practical  importance. 

---  For  low  noise  measurements  the  iterative  techniques* 
although  requiring  an  increase  in  computinq*  are  capable  of 
providing  better  tracking  for  maneuvering  targets. 

---  Tracking  with  only  one  buoy  is  acceptable  only  for 
non-maneuvering  targets. 

---  The  relative  position  of  the  buoys  is  one  of  the 
most  important  factors  which  influences  the  quality  of  the 


tracking 


VI.  CONCLUSIONS 


A.  SUMMARY  OF  RESULTS 

The  optimal  estimation  of  characteristics  and/or  parame- 
ters of  a certain  class  of  nonlinear*  dynamic*  stochastic 
systems  has  been  studied  in  a probabilistic  environment 
using  Bayes  formulation  concepts.  Approximate  solutions  and 
filtering  algorithms  were  generated  and*  amonq  them,  the 
known  Extended  Kalman  Filter  equations  and  higher  order 
filtering  equations  have  been  obtained  through  this  method. 

The  problem  of  tracking  submarine  targets  using  special 
passive  sonobuoys  was  modelled  and  with  this  model  extensive 
simulations  were  executed  to  allow  the  study  of  the  problem 
in  detail. 

Most  of  the  results  indicate  that  the  frequency  measure- 
ments have  minimal  effect  on  the  filtering  process.  The 
small  contribution  to  range  information  that  thev  do  pro- 
vide* when  associated  with  bearing  measurements*  can  nor- 
mally be  obtained  otherwise  by  judicious  placement  of  subse- 
quent buoys. 

The  utilization  of  other  types  of  measurements  such  as 
the  frequency  difference*  the  frequency  ratio  and  the  time 
delay  measurements*  proves  to  be  a great  help  in  improvinq 
computing  efficiency  by  eliminating  the  rest  frequency  as  a 


necessary  state*  and  thus  reducing  the  dimensionality  by 
one*  Also*  especially  with  the  time  delay  measurement*  a 
great  improvement  in  tracking  accuracy  is  possible. 

The  concept  of  partitioning  the  measurements  and"  pro- 
cessing them  separately*  even  if  they  occur  simultaneously* 
is  shown  to  bring  advantages  in  computing  efficiency  and 
also*  for  nonlinear  measurements*  in  tracking  accuracy. 
This  concept  is  of  great  practical  importance,  for,  in 
situations  where  the  measurements  are  sporadically  and  non* 
periodically  received  it  is  most  important  to  be  able  to 
process  the  measurements  as  they  naturally  occur. 

The  graphical  interpretation  of  the  action  of  Kalman 
filters*  developed  in  this  work*  provides  insight  into  the 
importance  of  each  variable  of  the  problem  in  the  filtering 
process.  The  direction  and  magnitude  of  the  correction  which 
is  applied  to  the  predicted  values  to  generate  the  new  esti- 
mate values  can  now  be  anticipated  as  a function  of  the 
measurement . 

Nonlinearity  errors  have  been  graphically  presented  and 
iterative  techniques*  including  the  known  Iterated  Extended 
Kalman  Filter  equations*  have  been  suggested  to  counteract 
their  disruptive  effect.  The  application  of  these  tech- 
niques to  the  tracking  problem  shows  that  improvement  in 
tracking  accuracy  is  possible. 

The  graphical  interpretation  also  indicates  the  very 
practical  conclusion  that  the  error  ellipses  tend  to  align 
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with  the  measurement  lines*  as  defined  in  Chapter  IV.  This 


provides  guidance  for  optimal  positioning  of  the  buoys  and 
the  types  of  measurements  to  process.  This  observation  is 
reinforced  by  the  simulation  results  of  Chapter  V. 


B.  SUGGESTIONS  FOR  FUTURE  RESEARCH 

The  concept  of  partitioning  the  measurements  allows  one 
to.  process  the  measurements  separately  and  opens  some  ques- 
tions ‘for  future  research*  such  as*  in  which  order  should 
nonlinear  measurements  be  processed  to  obtain  maximum  effi- 
c i ency? 

The  graphical  i nterpretat ion  of  the  filtering  process 
allows  anticipation  of  the  direction  and  magnitude  of  the 
corrections  which  are  then  applied  eo  predictions  to  gen- 
erate new  estimates.  This  situation  suggests  future  research 
in  the  determination  of  the  optimum  characteristics  of  the 
measurement  functions  which  in  turn  can  determine  the 
optimum  positions  and  characteristics  of  sensors. 

The  expansion  of  the  model  to  include  accelerations 
should  be  considered  . if  maneuvering  targets  have  to  be 
tracked  efficiently  and  the  increase  in  computing  power  can 
be  obtained.  Consideration  of  the  depth  of  the  tarqet  and 
the  inclusion  of  uncertainty  about  the  oosition  of  the  buoys 
are  additional  ways  to  extend  this  study. 


The  mode)  and  simulation  implemented  herein/  without 
benefit  of  classified  information#  is  very  flexible#  but 
many  idealizing  assumptions  have  been  made.  Ihus#  a natural 
extension  of  this  work  is  to  apply  the  algorithms  and  con- 
cepts developed  to  a real  problem  using  actual  sensor  data. 


APPENDIX  A : PROBABILITY  RELATIONS 


The  relation 

x = a + G . v (A.l) 

is  given  where  x and  v are  random  vectors  of  dimensions  n 
and  q*  r espec t i ve 1 y , a is  a n-d i mens i ona 1 deterministic  vec- 
tor and  G is  a Cnxq)  matrix  of  deterministic  coefficients. 

The  joint  density  Py(v.)  is  given  and  the  joint  density 
p x(  x_)  i s want ed . 


Case  ( i ) - If  n = q and  G 1 exists*  the  solution  is 
given  by  1133  * 

v_  - G 4 x^  - a^  = _f^ C_*  J (A. 2) 

and 


p C x ) 
x — 


P ( f C x ) ) 


(A. 3) 


Now  * 


and 


if  G,-i  exists# 
v = Grl  t x ’ - a * 1 = f (x1)  (A. a) 


From  Case  (i)»  the  solution  is 

p/*.') = ldet[4'^j  I • °v (i(-' } ) 


CA.5) 


P (x)  = p (x').dx  ,-dx  , ...dx  (A. o) 

x — x n+1  n+2  q 


Case  ( i i i ) *•  If  n > q the  following  development 


appl i es 


From  (A.l),  one  also  has  that 


= a’  ♦ G*  .v 

and  then*  if  G'  1 exists* 

v = G ,_1  [ x ' - a ' 1 = f ( x ’ ) 


From  Case  (i)  one  then  has 


p ix ,»x  ,...,x  ) 2 

x 1 2 q 


. D (f (x*  )) 
v — — 


Now  consider  the  variable  x Given  (A. 7), 

q+1 


p(x!x,x,...,x  ) = P(x!v,v,...*v  ) 

q+1  12  q q+1  12  q 


Out 


Vi=  Vi*  Vu'i  * V,2v2  * •••  * Wq 


= a . , + g _ .v 

q+1  -^q+l  — 


or 
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~ q+i  = Vi  * Vl  -- 


.ft*')  = 1 


(A. 9) 


From  (A. 9)  comes 


P l * .1  ! *,  # X-  » . . . , * ) - 6 ( x - 1 ) 

9+112  q q+1  'q+l 


P(x  ,x  ) = 6(x  - 1 ).p  (x  ,x  ,...,x  ) 

1 2 q+1  q+1  q+1  x'  1 2 q 


Doing  the  same  thing  for  x , x , one  finally 

q+^  n 


gets 


P (x) 
x — 


det  J_f 
3 xT 


p (f(x’)) 

v — — 


i=q+l 


(A.  10) 


For  the  SDecial  case  of  q = 1,  the  vectors  _x ' , a'  and 
9.'  become  scalars  and  the  sbove  equations  lead  to 

f Cx')  = f (x,  ) s (x  - a ) / g 
1 111 


= 1 / gi  ' 


1 = a +(x  -a)g/g 

r r 1 1 r 1 


and  Equation  (A. 10)  becomes 


P <x)  = ,5Cx  - 1 ) . 1 /g  . p C ( x - a )/g  ] 

x 11  lvll  l 


(A.  1 1 ) 


Another  similar  formula  can  be  obtained  for  this  sne* 
9 = 1 situation.  The  matrix  G in  (A.l)  is  now  a n- 
d i mens i ona 1 vector  and,  if  one  breaks  (A.l)  into  its  n 
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individual  scalar  equations  one  can  obtain 


pfx^  v) 


6 ( - at  - g^) 


Also  note  that 

ptx^iXj^v)  s ptxjv)  i P j 

and 


r\  ( ■ 


And  from  these 


re  1 at i ons  t 


Px/v(x!v)  =|| 
i=l 


p ( x ! v ) 


and  final  1 y » 


i=l 


p ( x ) 

X 


n 


i=l 


- g^l  .dv 


(A.  12) 


167 


APPENDIX  B : COMPONENTS  OF  COVARIANCE  MATRIX 


Fop  a single  measurement  of  the  form  z = H.j(  ♦ v»  the 
observation  matrix  becomes  the  vector 

Hs  lh  h •••  hi 

- 1 2 n 

where  n is  the  dimension  of  the  system. 


The  Kalman  gain  is  given  by 

G = P*HT[HP'HT  ♦ r2l  CB.l) 

where  P'  is  used  to  represent  the  prediction  error  covari- 
ance matrix,  which  is  symmetric,  and  r is  the  standard  devi- 
ation of  the  measurement  error. 

Since  H is  in  this  case  a vector,  the  product  HP'H 
will  be  a scalar. 


Let 


c = [HP'H1  ♦ r2) 


2,-1 


(B.2) 


The  vector  P'H  can  be  shown  to  be 


r A °LJ 

n 


P'H  = 


Z d’  h 
j-i  i 


Z d\  h 
j-1  nJ  j 


(6.3) 


and  thus  the  constant  c is  given  by 
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n n 2, 

C = l/  lE  I P h h ♦ r] 
i*l  j-1  ij  i j 


(B.4) 


From  (B.l)  one  now  sees  that  the  Kalman  gain  is  a vec- 
tor whose  comoonents  are  the  components  of  P'HT  multiplied 
by  the  constant  c,  i.e.  the  components  of  the  gain  vector 


• k si  lB'5) 


The  estimation  error  covariance  matrix  is  computed  from 


the  equation 


P r II  - G.H1  .P’  = P’  ♦ &P' 


(B.6) 


where 


AP*  = - GHP 1 = - cP’H  HP' 


(B.7) 


Let’s  define 


ai 


= o ’ h ♦ p ’ h ♦ 
il  1 i2  2 


,+p’h  (B.8) 

in  n 


Then » f rom  ( B . 3 ) » 


HP’  = (P’H  ) = [aj^  a2 


* an1 


(B . 9 ) 


and  from  (B.4), 


n o — 1 

c = llSl  * p 1 


(B.  10) 


From  (B.7)  the  increment  aP’  to  the  error  covariance 


matrix  is  then 
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APPENDIX  C : CHARACTERISTICS  OF  MEASUREMENT  FUNCTIONS 


in  Figure  C-l  some  curves  of  the  form  f ( jc ) s constant 
are  shown/  where  f(.x)  is  assumed  given.  The  vector  £a  and  a 
constant  value  d are  also  given.  The  vector  £b  is  defined 
by  the  equation 


xa  ♦ 


d h^ 

1*1  ’ III! 


(C-l) 


or 


b , 


x = 


xa  ♦ h.fd/.'h,'2] 


( C-2 ) 


where  h_  is  the  gradient  of  f(j^)  with  respect  to  x_  taken  at 


^ cl 

point  (x.  $ x_  )#  and  h/l^hj  is  a unit  vector  along  h . 


hi  s 4ft£a) 


I h | = (h2  ♦ h2)1/2 


h,  *7  f (xa) 
2 8X2  — 


lo  dete  mine  the  location  of  _x°  numerical  data  is 
required.  However/  it  is  possible  to  determine  qualitatively 
whether  x_b  ends  on  the  curve  f(jxa)  ♦ d/  or  on  one  side  or 
the  other  of  this  curve.  This  is  done  in  the  following 
paragraphs  for  three  special  functions. 


Case  (i)  - For  linear  f (jO  / j<D  always  ends  on  the  curve 
f(xj  =flxa)  ♦ d/  as  shown  below 
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r" 


r > 


S r 


Let  t ( jt ) = ♦ q*2 

I hen  h = p # h2  = q and  | h |2  = P2  ♦ q2 

h pom  Equation  (C-2), 

£b  = x_a  ♦ (pei  ♦ qe2).d/(p2  ♦ q2 ) 

so 

x b = x a t Pd/ (° 2 + q2 } 

* 2 = x a ♦ qd/ (d2  + q2) 

Now  t 

f(lb)  s p^  * q*2  5 pxi  + p2d/(p2  + q2)  + 

+ qxa  t q2  d/  (p2  + q2  ) = 

= p xa  ♦ q + d = f (xa)  ♦ d 

Case  ( i i ) - For  f(x_)  defined  as  below*  x_b  always 
on  top  of  f (_xa  ) ♦ df  as  in  the  Case  (i). 

f (j«  J = t C w1  - p)2  + (x2  - q )2  1 1/2 

Now  * 

hl  = (X1  " p)/m  h2  = (X2  " q)/m 

where 

m = l(xa  - p)2  ♦ (*2  - q)2l1/2=  f(*a> 

and 


ends 


.J 

j 
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From  (C-2)» 


«b  s x.a  ♦ h^.d  = _xa  ♦ l ( xa  - p)e^  ♦ C xa  " a)e_2^  *d^m 


so 


b _ a . , , _ a 

= *1  ♦ h^d  - x^ 

b _ a . . . «,  a 

*2  ' *2  * h2d  - x2 


♦ (i 


♦ ( 


- p)d/m 

- q)d/m 


Now  i 


fCxN  = l(xa  ♦ hxd  - P)2  + (>£  ♦ ^ d - q)2]^2^ 


s l(xa  * p)2  ♦ 2h  d(xa  • p)  + h2d2* 


♦ Cxa  - q)2  ♦ 2h  d(xa  - q)  ♦ h^d2!1^2 


Using  the  values  of  h^  and  h^» 


t(xb)  = l(xa 


d)2  ♦ (xa  - q)2  + d2  ♦ 


♦ 2dl(xa  - p)  /m  + ( xa  - q)  / m ] ^2 : 

1 2 


= (f(xa)2  + 2d.f(xa)  + d2]1/2= 


= [ C f ( x a)  + dJ2  J1/2= 


= f(xa)  ♦ d 


Case  C i i i ) * For  f(x)  as  defined  below  and  shown 
Figure  C-2»  x_b  always  ends  before  the  curve  f(^a)  ♦ d» 
[d j < v . 

f(x^)  = arctan((*2  - q)/(x^  - p)l 


i n 
for 
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Ihe  case  |dl  > it ✓ 2 is  not  of  interest  since  the  1 
and  fUJ  = f(xa)  ♦ d are  diverging.  For  |d|  <n/2» 


hx  = - { *2  " q)/n 


h = (xa  - p)/n 
2 1 


where 


and 


n s ( xa  • p)2  ♦ l*®  * 


| h | 2 = h2  ♦ h2  = 1/n 


(•  rom  Equation  (C-2)» 


x_b  = xa  t d.t-(xa  - qJe^  t (xa  - P)e2l 


so» 


b _ 

*1 

b _ 
x „ - 


a _ t -a 
x.  “lx 


1 

a 

*2 


o ) d 

p ) d 


Now  i 


t (xb  ) 


arctant(xb  - q)/(xb  - p)l  = 


or 


arctan  l(xa  - q)  ♦ ( xa  - p)dJ / 


[ ( xa  - o)  - ( x^  - q)dl 


t(xbJ  = arctanKs  ♦ d ) / ( 1 " sd)J 


where 


s = lx^  - q)/(*a  - P)  = tanlftx  )1 


i nes  h 
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tan[f(x_  ))  = (s  t d)/(l  ~ sd)  = 

= (tanlf(xa)i  ♦ dJ  / il  - d.tanlf  (xa)JJ 
= tanlfCxf  ) ♦ arctan(d)) 
or 

t (2Lb  ) - f (j«a  ) + arctan(d) 

Since  |d|  < tt/2,  then  arctan(d)  < d and  _xb  wil 
reach  the  line  f ( x a)  ♦ d. 


AHPtNOlX  D : ANALYTICAL  EVALUATION  OF  KALMAN  FILTERS 

After  the  design  of  a filter  an  analysis  phase  is 
necessary  to  verify  its  behavior  in  various  represent  at i ve 
situations. 

A "filtering  situation"  is  considered  here  to  be  com- 
pletely specified  by: 

---  the  true  parameters  and  initial  conditions  of  the 
system  which  generate  a unique  track  x(0),  x(l)»  ...»  x(t  ). 

---  the  true  parameters  of  the  sensors  and  their  meas- 
urement schedules. 

---  the  parameters  and  initial  conditions  assumed  by 
the  filter  for  the  system  and  for  the  sensors. 

lhe  approach  normally  used  to  determine  filter  behavior 
is  to  simulate  the  desired  situation  and  execute  hundreds  or 
thousands  of  Monte  Carlo  runs  and  compute  sample  statistics. 
The  most  useful  results  of  this  process  are: 

am(k)  - the  sample  mean  of  the  estimated  state  vector 
at  time  t^  obtained  from  m Monte  Carlo  runs. 

m m 

£ Ik)  = a (k)  - x.  C k ) - the  sample  mean  of  the 

estimated  error  vector  at  time  t,  . 

k 

m 

b Ik)  - the  sample  second  central  moment  of  the  state 

(or  error)  estimates  about  am  ( k ) » at  t,  . 

— k 


Ihe  objective  of  this  section  is  to  study  the  possibil- 
ity of  obtaining  values  of  £(k),  £(k)  and  b(k),  which  are 
approximately  the  values  of  ( k ) , em(k)  and  t)m ( k ) when  m is 
very  large*  without  the  use  of  the  time  consuming  Monte 
Carlo  operat ions. 

1 ---  Theoretical  Solution 

Unce  a "filtering  situation",  as  stated  above,  is 
defined,  the  operation  of  a Kalman  Filter  is  described  by 
the  recursive  equations  below  where  only  the  dependence  of 
each  term  on  the  previous  estimates  is  stressed. 

_S(k  + l S k)  s f_(£(k) ) (0.1) 

P(ktlik)  = <*,(_£  (k)).P(k).<t>T(^(k))  + 

♦ g(_x(k)  ) .Q(k)  .gT(xJk)  ) (0.2) 

\ 

(» (_x  ( k + 1 !k),P(k*l  }k))  = P(k  + i;k).HT(i(k  + l!k)). 

IH  ( i.(  k ♦ 1 J k)  ) .P(ktl  ! k)  .hT(i(k  + l | k)  ) + H(k))"1 

( D . 3 ) 

x_(  k + 1 ) = 5_(kfl|k)  ♦ G(5_(k  + 1 !k)  ,P(k  + l } k) ) . lh(x(k  + l ) ) + 


♦ v ( k + 1 ) 


h ( x ( k * 1 ! k ) ) ) 


(0.4) 


FU+1)  s [I  - G(£(ktlJk)#P(kel{k))*H(x_(k4'l{k))l,P(kflSk) 

(0.5) 

in  this  case  these  equations  describe  a deterministic 
process  were  it  not  for  the  random  measurement  noise  >/(k  + l)# 
the  only  random  input  to  the  "filter  dynamics"  since  a 
unique  track  is  considered.  The  initial  condition  of  the 
filter,  x^O)  and  P(0)  are  also  deterministically  given. 

tquations  (0.4)  and  (D.5)  can  be  rewritten  as 
x(k  + l)  s * (i/k  + 1 |k),P(kM  !k))  ♦ * ( x_(  k+  1 1 k ) , P ( k*  1 k ) ) . v ( k + 1) 
P(k+1)  = f (£(k+l |k),P(k+l Ik)) 

where  <p  r V and  r are  generally  nonlinear  matrix  functions 
of  the  variables  indicated. 

Considering  (0.1)  and  ( D . 2 ) one  can  then  write 

£(k+l)  = f_*  (x_(k),P(k))  + g*(i  (k)  ,P(k) ) .v_(ktl ) (0.6) 

P(ktl)  s h*(i(k),P(k))  CD. 7) 

where  f_*#  g*#  _h*  are  also  generally  nonlinear  matrix  func- 
tions. Prom  (0.6)  and  (D.7)  one  can  clearly  see  now  that# 
if  v is  a discrete  white  Gaussian  noise  process#  the  joint 
process  ( j<  # P 1 is  Markov. 

If  one  considers  a new  vector  y^  formed  with  the  n ele- 
ments of  x and  the  n.(nfl)/2  distinct  elements  of  the  (nxn) 
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symmetric  matrix  P#  y is  an  [n . ( n+ 3 ) /21 -d i mens i ona 1 process. 


Equations  (D.6)  and  (0.7)  can  then  be  combined  into 


y(ktl)  = A ( y ( k ) ) + B ( y ( k ) ) . v ( k + 1 ) 


(D.8) 


where  A is  a vector  function  and  B a R to  R m matrix  func- 


tion# m = n , (n  + 3 ) /2 . 


The  complete  behavior  of  the  filter  can  be  described  by 
the  probability  law  of  which  is  obtained  from 

p (j£(  k+1  ) ) = p (y  (k  + 1)  ijr(k)}.p(_y(k)).dy(k)  (D.9) 


and  p (_y(  k + 1 ) |^(  k ) ) is  obtained  from  p(v(k+l))  and  Equation 
(D.tt)  in  the  manner  shown  in  Appendix  A.  The  initial  value 
_y  ( 0 J is  deterministic  since  it  contains  the  initial  condi- 
tion of  the  filter  which  is  given. 


2 ---  The  Linear  Case 


tor  the  soecial  case  of  linear  dynamics  and  observa- 
tions# the  solution  can  be  obtained  in  a simpler  way.  In 
this  case  the  functions  <f>  » g and  H are  not  functions  of  the 
state  est i mat es-  and  so  neither  is  G. 


Equations  (D.l)  - (0.5)  can  be  grouped  now  into 


x ( k ♦ 1 ) = <t>  (x(k+l  !k),P(k+l#k))  ♦ ¥ ( P ( k ♦ 1 i k ) ) . v ( k + 1 ) 


P ( k ♦ l ) s r ( P ( k ♦ 1 ! k ) ) 


and  finally#  in  a recursive  way# 
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(D.  10) 


xCktl)  = f_*(  S^(  k ) ,P(k)  ) ♦ g*(P(k)).v_(k) 

P(k  + 1 ) = h*(P(k) ) CO.  1 1 ) 

From  (0.11)  it  is  clear  that  the  estimation  error 
covariance  matrix  is  now  a deterministic  process  and  can  be 
precomputed#  together  with  the  gains. 

Equation  CD. 10)  then  becomes 

x.Ck  + 1)  = f_*(i(k))  ♦ 2*(ktl).v(k  + l) 

* 

Prom  (D.l)  - (0.5)  the  values  of  _f  and  g*  can  be  found 
to  be#  for  this  linear  case# 

x_(  k ♦ 1 ) = II  ♦ G(k  + 1 ) . H ( k + 1 )1  .4,  (k)  ,_x  (k)  ♦ 

t G ( k + 1 ) .H(k+1 ) ,T(k+ 1 ) t G(k+l).v(k+l) 


or 


i(ktl)  = S(k  + l).x(k)  + F(k+1)  + G ( k ♦ 1 ) . v ( k ♦ 1 ) (0.12) 


where 

SCk)  = II  + G(k).H(k)l  A (k) 

P ( k ) = G(k).M(k).j[(k) 

Ihe  sought  after  behavior  of  the  filter  can  be 

described  by  the  moments  below#  which  agree  with  112) 


I(ktl)  = Elx(ktl))  = S(ktl).E  (x(k))  ♦ F(k+1)  * 


♦ G(k+l).Elv(k+l)l 


(0.13) 


e(  k + 1 ) = a(ktl)  - x(k+l)  = S ( k f 1 ) ,E 1 i ( k ) ) ♦ 


* (I  - G(kM).H(ktl)).x(ktl)  ♦ 


♦ G(k+l).E(v(k+l)J  = 


s SUM).e(k)  ♦ D(ktl)  .x(k+l)  - S(k+l).x_(k)  ♦ 


♦G(k*l).E(v(ktl)J 


(D.  1 a ) 


b ( k + 1 ) = S(ktl).b(k).Si(ktl)  + 


♦ G(kf 1 ) .E Lv(ktl ) .vT(k+l )) GT(k+l ) (0.15) 


3 ---  General  Case 


For  the  general  nonlinear  problem  one  returns  to  EQua- 
tion  (D.9;i  which  cannot  normally  be  solved  in  closed  form. 

Numerical  techniques  and  approximations  would  be  neces- 
sary that  may  use  more  computing  power  ana  present  worse 
results  than  the  Monte  Carlo  process  that  we  are  trying  to 


Une  way  to  simplify  the  problem  is  to  assume  that  y(k) 
has  an  approximate  Gaussian  distribution,  even  with  all  the 
non  1 i near i t i es  shown.  Linearizing  Equation  (D.8)  around  the 
mean  value  of  y_(  k ) would  give 
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1 


*(k+l)  s A(y(k))  ♦ A W(k) ) . [y  (k)  - y(k)J  + 


♦ B(y (k) ) .v(k+l ) 


(0.16) 


Taking  ^(0)  as  the  deterministically  known  initial 
value  y^O),  and  applying  the  expectation  and  variance  opera- 
tors to  Equation  (0.16),  one  gets  the  recursive  equations 
for  the  moments  of  y^  k ) : 


y.(ktl)  s A(y(k))  ♦ B(y(k)).Elv(k+l))  (0.17) 


Var  ly(k  + l)l  = I ~ A(y(k))l  .Varty(k))  . ( ~ A(y(k))lT  + 

3Z  - “ 3)1 

+ B(y(k) ) .Var tv(k  + l )]  .8T(y (k) ) (D.18) 

Ihe  moments  we  are  interested  in,  £,  e_  and  b are 

directly  obtained  from  subvectors  and  submatrices  of  the 
above  moments. 

The  primary  difficulty  with  this  approach,  however,  is 

in  obtaining  the  functions  A and  B.  This  can  be  seen  for  the 

very  simple  case  below,  as  an  example. 

Suppose  a scalar  system  with  a sinqle  observation  is 
described  by  the  equations 

x ( k+ 1 ) = x ^ ( k ) + w(k) 


z ( k ) 


x ( k ) + v ( k ) 


Equations  (0.1)  - (0.5)  give 


i(k  + l Ik) 

p(ktl ! k) 
l»  ( k+  1 ) = 


= x 2 ( k ) ; ^(k)  = <?x  ( k ) 

= 4i  2(k)  .pfk)  + q2  (k) 

ai2  (k)  .p(k)  ♦ q2 ( k ) 
ai2(k).p(k)  + q2 ( k ) + p2(k) 


and»  after  the  appropriate  stepsr 


x(ktl) 


p ( k 1 1 ) 


- 2 

xZ(k)  f 


^J2(k) .p(k)  + q2 (k) 
ai2(k).p(k)  + q2(k)  + r 2 ( k ) 
(ai2 (k)  .p(k)  + q 2( k ) ] .r2 (k) 
ax2(k)  .p(k)  + q2 (k)  + r2 (k) 


.Ix(k)  - x(k)J 


ai2(k).p(k)  + q2 ( k ) 
ai2(k).p(k)  ♦ q 2 ( k ) + r2 ( k ) 


. v ( k + 1 ) 


Ihe  mean  value  given  by  Equation  (D.17)  can  be  obtained 
in  a simpler  way » however.  It  is  the  result  obtained  by 
running  the  filter  in  the  simulated  environment  but  making 
the  measurement  noise  assume  its  mean  value  E(v(k))»  nor- 
mally zero . 


If  one  uses  more  terms  in  the  expansion  of  Equation 
(0.8)#  better  results  will  be  obtained  at  the  expense  of 
increased  computing  effort  and  time. 
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P or  this  general  case  it  is  suggested  that  these  equa- 
tions be  applied  to  a simple  scalar  or  two-dimensional  non- 
linear problem  and  that  the  computing  power  required  to  find 
the  sought  after  moments  be  compared  with  that  required  to 
run  a Monte  Carlo  process  yielding  the  same  level  of  accu- 
racy in  the  results. 


APPENOIX  E : COMPUTER  PROGRAM 


Ihe  basic  structure  of  the  computer  program  written  for 
the  simulation  of  the  tracking  problem*  and  the  evaluation 
of  models  and  filtering  algorithms*  is  shown  in  Figure  E-l. 
The  solid  lines  represent  the  normal  flow  of  control  within 
the  program*  the  dashed  lines  show  the  extra  pathes  that 
become  available  whenever  an  interruption  is  requested  by 
the  user. 

After  a brief  introduction  to  the  program  a table 
called  MtNU  is  presented  to  the  user*  as  shown  in  Figure  E- 
2.  If  the  user  chooses  actions  number  1 or  5 the  result  is 
clear. 

Choice  number  3 provides  access  to  all  the  problem 
variables.  Since  the  number  of  variables  that  characterize 
each  simulation  is  very  large*  a simple  way  to  deal  with 
these  variables  had  to  be  devised.  As  it  is  implemented* 
the  program  always  "remember"  the  values  given  to  the  vari- 
ables at  the  last  time  the  program  was  used*  so  that  only 
the  variables  to  have  their  values  modified  have  to  be 
addressed.  The  change  of  value  of  a variable  is  simply  maae 
by  selecting  the  page  where  it  appears*  typing  the  letter 
associated  with  thevariable  and  the  new  desired  value.  The 
program  immediately  responds  by  presenting  back  the  entered 
value.  A more  detailed  diagram  of  the  actions  resulting 


KFNti 


f 


* 


1 


Press  the  rnjTCer  corresoondi  ng  to  the  action 
desired  and  c / r . 

\ - PRFSE  NT  TUTORIAL  INFORMATION 

2 - MODIFY  PROGRAM  FLAGS 

3 - FORMULA  IF  OH  MOD  T F Y THE  PROBLEM 

a - STAHT  thf  SOIUTION  OR  THF  PROBLEM 
5 - END  THF  PROGRAM  AND  RXTI 


Fiaure  E-2  I I Re  MENU 
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from  choice  number  3 is  presented  in  Figure  E-3.  A typical 
page  of  variables  is  shown  in  Figure  E-9 . 


* 

r 


Choice  number  4 simulates  the  tracking  problem  accord- 
ing to  the  values  defined  in  the  previous  step.  A block 
diagram  of  the  actions  involved  in  the  simulation  is 
presented  in  Figures  E-5 , E-6  and  E-7. 

Initially  all  the  important  problem  variables  are 
printed  with  their  current  values  in  a form  as  shown  in  Fig- 
ures fc-B,  E-9  and  E-10.  During  the  first  run  the  scheduling 
of  all  the  events  involved  in  the  simulation  is  also  printed 
for  future  reference*  as  shown  in  Figure  E-ll. 

with  choice  number  2 from  the  MENU  a new  table  is 
presented  as  shown  in  Figure  E-12.  This  table  in  also 
presented  to  the  user  whenever  he  requests  an  interruption, 
at  any  time  or  point  within  the  program. 

The  extra  choices  that  now  become  available  are  almost 
self-explanatory  but  it  should  be  added  that  with  choice 
number  7 the  simulation  is  ended  and  the  partial  statistics 
are  computed  and  written  in  the  appropriate  output  files; 
choice  number  6 allows  the  modification  of  any  variable  in 
the  middle  of  the  simulation,  in  the  same  way  as  explained 
before*  choice  number  0 transfer  control  back  to  MENU  or  to 
the  point  of  i nterrupt i on. 

A listing  of  the  program,  in  C language,  is  available 
at  the  Electrical  Engineering  Department,  Naval  Postgraduate 


190 


To  change  the  value  of  any  variable*  type  the  indicated 
letter,  at  least  one  space,  the  desired  value  and  c/r. 

To  see  the  other  variables  of  the  problem,  type  1 c/r. 

To  qo  back  to  menu,  tyoe  0 c/r, 

( I meter  = t.09<l  yards;  1 mc-ter/sec  - 1 «0<*  knots  ) 


TARGFT  INITIAL  PARAMFTFRS 


initial  * - position 

: 0.0 

km 

a 

y - position 

: 0.0 

km 

b 

speed 

: 7.0 

m/sec 

c 

head i no 

: 30S.0 

deq 

d 

f reauency 

: S00.A 

hertz 

e 

standard  deviation  of 

f ore i no 

func tions 

i 

speed/sec 

! 0.0 

m/sec/s®c 

f 

Headi no/sec 

: o.o 

deq/sec 

q 

f reauenc v/sec 

t 0.0 

hert  z/sec 

h 

Figure  E-A 


Typical 


oaoe 


of  variables 


from  MFNU 


Print  Val ue  of 
AM  Variables 


run  = t 


Initialize  Variables 
and  Random  Number  Gen, 
for  a New  Run 


time  = 1 


run  = 1 ? 


Is  taroet  to  be  updated 
everv  second? 


l ? 


Are  There  Anv 
Other  Fvents 
Occurrino  at 
This  Time? 


Print  Table 
of  Events 


! Update 
■ Target 
; Parameters 


/?> 


Figure  E-b  : Simulation  T 


Is  time  > durat i on 
of  a run?  - — — 


Is  run  > number 
o f runs?  — " 


Compute 

Statistics 


Send  Kesu 1 1 s 
To  OutDUt 

Piles 


f o 

MENU 


Fioure  E-b  : Simulation  II 


Select  Action 
From  Table 


Chanae  Taroet 
Parameters  to 
Simulate 
Maneuver 


Store  Values 
of  Parameters 
and  Fstimates 
for  Statistical 
Purposes 


Is  taroet  to  be 
updated  before 
measurements? — ■ 


t'odate 

Target 

Parameters 


Simulate 
Sel ected 
Measurement 


Determine 

Predicted 

State 


L i near i ze 
About 
Predicted 
State 


Process 
Measurement 
To  Obtain 
Fst i mate 


Is  the  requested 
number  of  iterations 
sat i s f i ed? 


t i n e a r i z e 
About  Partial 
Fst imate 


Fioure  E-7  : Simulation  TII 


TARGET  INITIAL  PAHAMFTFRS 


initial  * - Dosition 

0.0 

km 

a 

v “ oos i t i on 

0.0 

km 

b 

speed 

7.0 

m/sec 

c 

s 

headi no 

305.0 

deq 

d 

1 

f requenc v 

500. A 

hertz 

e 

standard  deviation  of 

Pore i na 

f unct i ons 

speed/sec 

0.0 

m/sec/sec 

f 

headi na/sec 

0.0 

deq/sec 

o 

f requenc v/sec 

0.0 

hert  z / sec 

h 

GENERAL  DIRECTIONS 


number  of  PIFAR  buoys  : 5 

buoys 

a 

number  of  L0FAR  buoys 

: 0 

buovs 

b 

number  of  maneuvers 

2 

man 

c 

number  of  runs 

?0 

runs 

d 

durat ion  of  runs 

2000 

Sec 

e 

number  of  key  ooints 

20 

points 

f 

seed 

13339 

o 

average  sound  velocity 

1500.0 

m/sec 

h 

range  of  the  buovs 

8.0 

km 

i 

FILTER  INITIAL  PARAVFTFRS 


initial  * • Position  ; -0.5  k"  a 
V - Position  .!  . 0 km  b 

«o^eH  : 6.0  m/s*»c  c 
••eaciino  : MO.rt  d«»r;  d 
frequency  : 500.5  hortz  e 

standard  deviation  of  forcino  ♦unctions 

spreo/see  : O.OOS  m/sec/sec  f 
headino/sec  : 0.050  oeq/sec  o 
♦ requencv/sec  : 0.010  hertz/sec  h 


Fiaure  E-b  : listing  of  problem  variables  T 
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initial  covariance  matrix 


* DOS 

V POS 

speed 

heaoo 

f reo 

X POS 

UO  a 

0.3  b 

0.0  c 

0.0  d 

0.0  e 

y pos 

* 

1.0  q 

0.0  h 

0.0  i 

0.0  j 

speed 

* 

* 

2.0  m 

0.0  n 

0.0  o 

Heada 

* 

* 

* 

100.0  s 

0.0  t 

f req 

* 

* 

* 

* 

1.0  y 

( km , 

m/sec  # 

decrees  and 

hertz 

c rossmu 1 t i o 1 i ed) 

MEASUREMENTS  SCHFOI'LE  I 


huoy 

t yoe 

start 

period 

mean 

stdev 

1 

B 

100 

1 00 

0.0 

5.0 

abcaef 

1 

F 

100 

100 

0.0 

0.04 

chi j k 1 

? 

6 • 

500 

100 

0.0 

5.0 

mnnpqr 

2 

F 

500 

1 00 

0.0 

0.04 

stuvwx 

(start  and  period  in  seconds) 

(mean  and  std  dev  in  deqrees»  hertz  or  seconds) 
(period  must  oe  zero  if  ’reasurement  is  not  used) 


measurements  SCHEDULE  II 


buoy 

type 

start 

period 

mean 

stdev 

1? 

T 

500 

100 

0.0 

0.01 

abcdef 

13 

T 

500 

1 00 

0.0 

0.01 

oh i j k 1 

23 

T 

500 

1 00 

0.0 

0.01 

mnopor 

3 

B 

500 

1 00 

0.0 

5.0 

s t li  v w X 

(star t 

and 

pe  r i od 

in  seconds) 

(mean  and  std  dev  in  deqrees»  hertz  or  seconds) 
(period  must  oe  zero  if  measurement  is  not  used) 


Eioure  E - 9 * Listinq  of  problem  variables  II 
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TARGET  MANEUVERS 


1st  maneuver  : time 

T 

f>00 

sec 

a 

s • 

vdot 

- 

0.0 

m/sec /m i n 

b 

• 

hdot 

- 

-9.0 

deg/mi n 

c 

f dot 

r 

0.0 

hertz/mi n 

d 

TARGET  MANEUVERS 


2nd  maneuver  S time 

s 

1200 

sec 

a 

vdot 

r 

0.0 

m/sec/mi n 

b 

hdot 

5 

o.o 

deg/m i n 

c 

f dot 

r 

0.0 

hertz/m  in 

d 

RUQYS  PARAMFTFRS 

First  buoy  s type  - OIFAR  a 

x - position  - 0.0  km  b 

y - position  = -3.0  km  c. 

If  DIFAR: 

bearing  error  means  0.0  deg  d 

std  oev  = S.O  oeg  e 

a 

h 
c 


If  DTFAR: 

bearing  error  means  0.0  deg 

std  oev  = S.o  o »g 

Fioure  E-10  t listing  of  problem  variables  III 
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BUOYS  PARAMFTFRS 

2nd  buoy  : tvpe  = DTFAR 

v - position  - 4.0  km 

y - position  s -3.0  km 


i a 


SEQUENCE  OF  EVENTS  FRO*'  601  TO  900  SECONDS 


: . s 


v ' 


700  sec  “ bearino  measurement  by  buoy  number  1 
700  sec  ~ frequency  measurement  by  buoy  number  1 
700  sec  “ bearino  measurement  by  buov  number  2 
700  sec  “ frequency  measurement  by  buov  number  2 
700  sec  - time  delay  measurement  bv  buoys  1 and  2 
700  sec  • Monte  Carlo  coint  number  7 
#00  sec  - bearino  measurement  by  buov  number  1 
000  sec  “ frequency  measurement  by  buov  number  1 
000  sec  • bearino  measurement  by  buov  number  2 
000  sec  • frequency  measurement  by  buov  number  2 
000  sec  “ time  delay  measurement  by  buoys  1 and  2 
000  sec  - Monte  Carlo  noint  number  7 
900  sec  • bearino  measurement  by  buov  number  1 
900  sec  - frequency  measurement  by  buov  number  1 
900  sec  * bearino  measurement  by  buov  number  2 
°00  sec  “ f reauency  measurement  by  buov  number  2 
900  sec  - time  delay  measurement  by  buoys  1 and  ? 
900  sec  * Monte  Carlo  coint  number  7 

Fjoure  E-ll  : Table  of  events. 


— ■ — -- 
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Press  the  number  corresoondinq  to  the  option  or 
action  desired  and  c/r. 

To  continue  or  return  oress  0 c/r. 


1 - update  tercet  every  second 


? - update  taroet  only  before  neasurements 


3  “ start  orintinq  on  line  printer 


4  “ stoo  orintino 


5  ~ nut  results  of  next  run  on  output  files 


6  • modify  parameters  of  the  problem 


7  - terminate  the  problem 


.4 

s 

Fioure  E-l?  : Interruption  table. 
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